Katherine Analysis
- Housekeeping
- Introducing the Data
- Reviewing of hypotheses:
- Field Data Analysis
- Soil Characteristics
- Histology data Analysis
- Histology Data exploration
- Molecular data Analysis
Housekeeping
Preamble
By Katherine Duchesneau
We are working with the roots Lesion and Mycorrhiza dataset for now.
R details
This section contains all the necessary details about the R version and package version used to run this Notebook.
- File creation date: 2018-02-12
- R version 3.4.1 (2017-06-30)
tidyversepackage version: 1.2.1lme4package version 1.1.13MASSpackage version 7.3.47MuMInpackage version 1.15.6emmeanspackage version 0.9.1bootpackage version 1.3.20brmspackage version 2.0.1loopackage version 1.1.0fitdistrpluspackage version 1.0.9mefapackage version 3.2.7library(devtools)andinstall_github("vqv/ggbiplot")package version 0.55veganpackage version 2.4.3glmmADMBpackage version 0.8.3.3carpackage version 2.1.6factoextrapackage version 1.0.5.999pvclustpackage version 2.0.0
Load Packages and Functions
library(tidyverse)
library(lme4)
library(MASS)
library(MuMIn)
library(emmeans)
library(boot)
library(brms)
library(loo)
library(fitdistrplus)
library(mefa)
library(ggbiplot)
library(pvclust)
library(vegan)
library(ggbiplot)
library(factoextra)
library(glmmADMB)
library(car)Load a custome theme for clean, readable graphs:
Overdispersion function, tests for overdispersion in GLMM:
Introducing the Data
Load/ Modify data: Co-occurence
## [1] 99.55752## [1] 100## indiv2 None.Path None.Myc Herbivory
## 10IAH1 : 100 Length:22500 Length:22500 Min. :0.00000
## 10IAH2 : 100 Class :character Class :character 1st Qu.:0.00000
## 10ISC1 : 100 Mode :character Mode :character Median :0.00000
## 10ISC2 : 100 Mean :0.03658
## 10ISC3 : 100 3rd Qu.:0.00000
## 10OAH1 : 100 Max. :1.00000
## (Other):21900
## None.Herb Lesion Mycorrhiza
## Min. :0.0000 Min. :0.000 Min. :0.0000
## 1st Qu.:1.0000 1st Qu.:0.000 1st Qu.:0.0000
## Median :1.0000 Median :1.000 Median :1.0000
## Mean :0.9631 Mean :0.622 Mean :0.5732
## 3rd Qu.:1.0000 3rd Qu.:1.000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.000 Max. :1.0000
## ## indiv2 None.Path None.Myc Herbivory
## Length:22500 Min. :0.000 Min. :0.0000 Min. :0.00000
## Class :character 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:0.00000
## Mode :character Median :0.000 Median :0.0000 Median :0.00000
## Mean :0.378 Mean :0.4268 Mean :0.03658
## 3rd Qu.:1.000 3rd Qu.:1.0000 3rd Qu.:0.00000
## Max. :1.000 Max. :1.0000 Max. :1.00000
##
## None.Herb Lesion Mycorrhiza pop
## Min. :0.0000 Min. :0.000 Min. :0.0000 5 :2900
## 1st Qu.:1.0000 1st Qu.:0.000 1st Qu.:0.0000 7 :2600
## Median :1.0000 Median :1.000 Median :1.0000 8 :2600
## Mean :0.9634 Mean :0.622 Mean :0.5732 3 :2200
## 3rd Qu.:1.0000 3rd Qu.:1.000 3rd Qu.:1.0000 4 :2200
## Max. :1.0000 Max. :1.000 Max. :1.0000 1 :2000
## (Other):8000
## location species
## I:10600 AH:3100
## O:11900 CL:3900
## GA:4900
## GR:1300
## MR:3100
## SC:6200
## ## indiv2 None.Path None.Myc Herbivory None.Herb Lesion
## Length:22500 0:13995 0:12897 0:21677 0: 823 0: 8505
## Class :character 1: 8505 1: 9603 1: 823 1:21677 1:13995
## Mode :character
##
##
##
##
## Mycorrhiza pop location species
## 0: 9603 5 :2900 I:10600 AH:3100
## 1:12897 7 :2600 O:11900 CL:3900
## 8 :2600 GA:4900
## 3 :2200 GR:1300
## 4 :2200 MR:3100
## 1 :2000 SC:6200
## (Other):8000## [1] 225Load/ Modify data: Totals
## [1] 226## [1] 226## indiv Decay Pathogen Hyphae
## 10IAH1 : 1 Min. : 0.000 Min. : 0.000 Min. : 2.00
## 10IAH2 : 1 1st Qu.: 1.000 1st Qu.: 2.000 1st Qu.:34.00
## 10ISC1 : 1 Median : 4.000 Median : 4.000 Median :48.00
## 10ISC2 : 1 Mean : 7.894 Mean : 7.168 Mean :47.23
## 10ISC3 : 1 3rd Qu.: 9.000 3rd Qu.:10.000 3rd Qu.:63.00
## 10OAH1 : 1 Max. :60.000 Max. :86.000 Max. :96.00
## (Other):220
## None.Path Arbuscules Vesicles M_Hyphae
## Min. : 0.00 Min. : 0.000 Min. : 0.000 Min. : 0.00
## 1st Qu.:20.00 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.:30.25
## Median :33.00 Median : 0.000 Median : 6.000 Median :46.00
## Mean :37.71 Mean : 2.363 Mean : 9.956 Mean :45.16
## 3rd Qu.:53.75 3rd Qu.: 3.000 3rd Qu.:14.000 3rd Qu.:61.00
## Max. :97.00 Max. :21.000 Max. :73.000 Max. :94.00
##
## None.Myc Herbivory None.Herb pop location
## Min. : 0.00 Min. : 0.000 Min. : 62.00 5 :29 I:106
## 1st Qu.: 18.00 1st Qu.: 0.000 1st Qu.: 96.00 8 :27 O:120
## Median : 38.50 Median : 1.000 Median : 99.00 7 :26
## Mean : 42.52 Mean : 3.624 Mean : 96.38 3 :22
## 3rd Qu.: 66.00 3rd Qu.: 4.000 3rd Qu.:100.00 4 :22
## Max. :100.00 Max. :38.000 Max. :100.00 1 :20
## (Other):80
## species
## AH:31
## CL:39
## GA:50
## GR:13
## MR:31
## SC:62
## Field survey data
FloristicSurvey <-read.csv("CSV/May_2017_FloristicSurvey_summer2016.csv", sep="\t")
PlantLength <- read.csv("CSV/May_2017_PlantLength_summer2016.csv", sep="\t")
Soil_characteristics <- read.csv("CSV/Soil_characteristics.txt", sep="\t")
#Restructuring
head(PlantLength)str(PlantLength)## 'data.frame': 231 obs. of 4 variables:
## $ Pop : int 14 14 14 14 14 14 14 14 14 14 ...
## $ Location : Factor w/ 2 levels "in","out": 2 2 2 2 2 1 1 1 1 1 ...
## $ Species : Factor w/ 7 levels "Acer_saccharum",..: 1 4 3 5 7 1 4 3 5 7 ...
## $ Length_plant: num 14 87 20 12 67 11 65 51 20 100 ...colnames(PlantLength)[1] <- "Population"
PlantLength$Population<-as.factor(PlantLength$Population)
PlantLength$Location<-as.factor(PlantLength$Location)
FloristicSurvey$Population<-as.factor(FloristicSurvey$Population) # Make the population a factor and not an interger
FloristicSurvey[is.na(FloristicSurvey)] <- 0
sum(row.has.na <- apply(FloristicSurvey, 1, function(x){any(is.na(x))}))## [1] 0Root size
Root_size <-read.csv("CSV/Root_size.csv")
Root_size$pop<-gsub("^([0-9]+)[I,O][A-Z]+.*","\\1",Root_size$indiv)
Root_size$location<-gsub("^[0-9]+([I,O])[A-Z]+.*","\\1",Root_size$indiv)
Root_size$species<-gsub("^[0-9]+[I,O]([A-Z]+).*","\\1",Root_size$indiv)
Root_size$location<-as.factor(Root_size$location)
Root_size$species<-as.factor(Root_size$species)
Root_size$pop<-as.factor(Root_size$pop)Variable descriptions
The variables are:
Co-occurence:
indiv2: The unique code given to a sample. In this dataset the unique sample code repeats 100 times for each cross where an observation was taken on the root sample.
None.Path: A binary representation of whether a sign of pathogen activity was recorded (0) or not (1).
Lesion: The counterpart to the previous variable (None.Path). A binary representation of whether a sign of pathogen activity was recorded (1) or not (0).
None.Myc: A binary representation of whether a sign of myccorhizal activity was recorded (0) or not (1).
Mycorrhiza: The counterpart to the previous variable (None.Myc). A binary representation of whether a sign of mycorrhizal activity was recorded (1) or not (0).
None.Herb: A binary representation of whether herbivory was recorded (0) or not (1).
Herbivory: The counterpart to the previous variable (None.Herb). A binary representation of whether Herbivoryn was recorded (1) or not (0).
Population (pop): The coding number representing the population at which the sample was collected.
location: A code representing the whether the sample was collected inside a Alliaria petiolata population (I) or whther it was collected at least 7 m outside of he furthest individual in the A. petiolata population (O).
species: The particular species to which the sample belongs.
Cross: The particular cross number where the data was recorded on the individual sample.
Total scores:
indiv: The unique code given to a sample.
Decay: The total number of crosses where decay was recorded an individual sample when doing a total of 100 crosses/ sample.
Pathogen: The total number of crosses where pathogen was recorded on an individual sample when doing a total of 100 crosses/ sample.
Hyphae: The total number of crosses where non-myccorhizal hyphae was recorded on an individual sample when doing a total of 100 crosses/ sample.
None.Path: The total number of crosses where no signs of pathogen activities were recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Decay, Pathogen, Hyphae, and None.Path must come to 100 per individual to account for all 100 crosses.
Arbuscule: The total number of crosses where an arbuscule was recorded on an individual sample when doing a total of 100 crosses/ sample.
Vesicules: The total number of crosses where a vesicule was recorded on an individual sample when doing a total of 100 crosses/ sample.
M_Hyphae: The total number of crosses where myccorhizal hyphae was recorded on an individual sample when doing a total of 100 crosses/ sample.
None.Myc: The total number of crosses where no signs of mycorrhizal activities were recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Arbuscules, Vesicules, M_Hyphae, and None.Myc must come to 100 per individual to account for all 100 crosses.
Herbivory: The total number of crosses where herbivory was recorded on an individual sample when doing a total of 100 crosses/ sample.
None.Herb: The total number of crosses where no signs of herbivory was recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Herbivory, and None.Herb must come to 100 per individual to account for all 100 crosses.
Population (pop): The coding number representing the population at which the sample was collected.
location: A code representing the whether the sample was collected inside a Alliaria petiolata population (I) or whther it was collected at least 7 m outside of he furthest individual in the A. petiolata population (O).
species: The particular species to which the sample belongs.
Reviewing of hypotheses:
Question 1: Is there evidence that Mycorrhizal colonization provides protection against root pathogen damage?
Question 2: Do these effects vary among species?
Hypothesis 1
Mycorrhizal colonization protects roots from lesions and herbivory. 
Hypothesis 2
A. petiolata invasion reduces mycorrhizal colonization and increases pathogen colonization and herbivory 
Hypothesis 3
Cycorrhizal colonization is reduced in the presence of A. petiolata 
Hypothesis 4
Mycorrhizal protection is completely lost in the presence of A. petiolata 
Molecular data Hypothesis
Hypothesis 1
A. petiolata invasion changes the diversity of mycorrhizal and other fungal species present in the soil.
I will find that the molecular data from the soil extraction supports this claim by having a reduced biodiversity and abundance of mycorrhiza in samples taken inside A. petiolata populations.
Hypothesis 2
There will be a negative correlation between mycorrhizal and fungal pathogen diversity.
Field Data Analysis
Plant length and GM presence
row.has.na <- apply(PlantLength, 1, function(x){any(is.na(x))})
sum(row.has.na) ## [1] 2PlantLength.filtered <- PlantLength[!row.has.na,]# Removed all the rows with NAs
PlantLength.filtered$Pop <- as.factor(PlantLength.filtered$Pop)
car::qqp(PlantLength$Length_plant, "norm")
shapiro.test(PlantLength$Length_plant)##
## Shapiro-Wilk normality test
##
## data: PlantLength$Length_plant
## W = 0.89321, p-value = 1.14e-11mod1<-lm(Length_plant~Location*Species, data=PlantLength.filtered)
bc<-boxcox(mod1, lambda = seq(-2, 2, 0.1))
lambda <- bc$x[which.max(bc$y)]
PlantLength.filtered <- cbind(PlantLength.filtered, ((PlantLength.filtered$Length_plant^lambda)-1)/lambda)
names(PlantLength.filtered)[length(PlantLength.filtered)] <- "Yprime"
mod1BC <- lm(Yprime~Location*Species, data=PlantLength.filtered)
par(mfrow = c(2, 2))
plot(mod1BC)
# This is basically a log transformation.
# Multicolinearity?
car::vif(mod1BC)## GVIF Df GVIF^(1/(2*Df))
## Location 5.43980 1 2.332338
## Species 94.60756 6 1.461035
## Location:Species 306.05628 6 1.611202sqrt(car::vif(mod1BC)) > 2## GVIF Df GVIF^(1/(2*Df))
## Location TRUE FALSE FALSE
## Species TRUE TRUE FALSE
## Location:Species TRUE TRUE FALSEmod1<-lmer(log(Length_plant)~Location*Species+(1|Pop), data=PlantLength.filtered)
mod2<-lmer(log(Length_plant)~Location+(1|Pop), data=PlantLength.filtered)
mod3<-lmer(log(Length_plant)~1+(1|Pop), data=PlantLength.filtered)
AICres <- model.sel(mod1,mod2,mod3)
r.squaredLR(mod1)## [1] 0.7640355
## attr(,"adj.r.squared")
## [1] 0.8900227summary(mod1)## Linear mixed model fit by REML ['lmerMod']
## Formula: log(Length_plant) ~ Location * Species + (1 | Pop)
## Data: PlantLength.filtered
##
## REML criterion at convergence: 161.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.94568 -0.62634 -0.05365 0.56828 2.80893
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pop (Intercept) 0.006596 0.08122
## Residual 0.100384 0.31683
## Number of obs: 229, groups: Pop, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.34167 0.07565 30.955
## Locationout -0.03789 0.09778 -0.388
## SpeciesAnemone_hepatica 0.23682 0.12794 1.851
## SpeciesCircaea_lutetiana 1.08300 0.10771 10.055
## SpeciesGallium_aparine 1.46347 0.10204 14.342
## SpeciesGeranium_robertinum 0.80086 0.10771 7.435
## SpeciesMaianthemum_racemosum 0.55442 0.12295 4.509
## SpeciesSolidago_canadensis 1.63655 0.09778 16.738
## Locationout:SpeciesAnemone_hepatica 0.21256 0.17101 1.243
## Locationout:SpeciesCircaea_lutetiana -0.24113 0.14797 -1.630
## Locationout:SpeciesGallium_aparine -0.20452 0.14392 -1.421
## Locationout:SpeciesGeranium_robertinum 0.07377 0.14799 0.498
## Locationout:SpeciesMaianthemum_racemosum 0.12277 0.16743 0.733
## Locationout:SpeciesSolidago_canadensis -0.11605 0.13828 -0.839##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itconfint(mod1)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 0.01963921 0.16258131
## .sigma 0.28096327 0.33850558
## (Intercept) 2.19817644 2.48515490
## Locationout -0.22469501 0.14891543
## SpeciesAnemone_hepatica -0.01319870 0.47949790
## SpeciesCircaea_lutetiana 0.87828687 1.29033635
## SpeciesGallium_aparine 1.26920608 1.65925061
## SpeciesGeranium_robertinum 0.59591108 1.00763645
## SpeciesMaianthemum_racemosum 0.31810708 0.78817334
## SpeciesSolidago_canadensis 1.44974225 1.82335269
## Locationout:SpeciesAnemone_hepatica -0.11265649 0.54145529
## Locationout:SpeciesCircaea_lutetiana -0.52495034 0.04066983
## Locationout:SpeciesGallium_aparine -0.47948584 0.07045383
## Locationout:SpeciesGeranium_robertinum -0.20959317 0.35592002
## Locationout:SpeciesMaianthemum_racemosum -0.19718336 0.44252670
## Locationout:SpeciesSolidago_canadensis -0.38023461 0.14813034Floristic composition and GM presence
Linear discriminate analysis
# Discriminant function analysis is used (general LDA format) because i have a cat. dependent and a binary independent. (it needs binary or continuous!)
# you might have to bin covers because it only works with categorical, so 0-5% none, 6-40 low etc.
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("None", "Med_Low", "Med_High", "Overtaken"))
table(FloristicSurvey$GM_Coverage_category)##
## None Med_Low Med_High Overtaken
## 82 24 58 34# [0-20] None
# ]20-50] Medium low
# ]50-80] Medium high
# ]80-100] Overtaken#LDA
FloristicSurveysmall<-data.frame(FloristicSurvey[,4:56])
# Linear discriminant function analysis
FloristicSurveysmall <- FloristicSurveysmall[,-c(42,45)]
FloristicSurveyLDA <- lda(GM_Coverage_category ~., data=FloristicSurveysmall)
# Extract scaling vectors
scalvec<-data.frame(FloristicSurveyLDA$scaling)
# Extract predictions
FloristicSurveyLDAval <- data.frame(predict(FloristicSurveyLDA)$x)
# Plot results
ggplot(data=FloristicSurveyLDAval,aes(x=LD1, y=LD2, group=FloristicSurveysmall$GM_Coverage_category))+
stat_ellipse(geom="polygon",aes(colour=FloristicSurveysmall$GM_Coverage_category),fill=NA,size=1.2,alpha=0.3)+
stat_ellipse(geom="polygon",aes(fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=1.2,alpha=0.3)+
geom_point(aes(shape=FloristicSurveysmall$GM_Coverage_category,fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=I(4),alpha=I(0.8))+
xlab("LD Axis 1")+ylab("LD Axis 2")+theme_simple() +
theme(legend.title=element_blank())
#test significance of axis
anova(lm(FloristicSurveyLDAval$LD1~FloristicSurveysmall$GM_Coverage_category))anova(lm(FloristicSurveyLDAval$LD2~FloristicSurveysmall$GM_Coverage_category))#LDA WITHOUT A. PETIOLATA
FloristicSurveysmall<-data.frame(FloristicSurvey[,4:56])
# Linear discriminant function analysis
FloristicSurveysmall[10]<-NULL
FloristicSurveysmall[c(41,44)] <- NULL
FloristicSurveyLDA<-lda(GM_Coverage_category ~ ., data=FloristicSurveysmall)
# Extract scaling vectors
scalvec<-data.frame(FloristicSurveyLDA$scaling)
# Extract predictions
FloristicSurveyLDAval <- data.frame(predict(FloristicSurveyLDA)$x)
# Plot results
ggplot(data=FloristicSurveyLDAval,aes(x=LD1,y=LD2,group=FloristicSurveysmall$GM_Coverage_category))+
stat_ellipse(geom="polygon",aes(colour=FloristicSurveysmall$GM_Coverage_category),fill=NA,size=1.2,alpha=0.3)+
stat_ellipse(geom="polygon",aes(fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=1.2,alpha=0.3)+
geom_point(aes(shape=FloristicSurveysmall$GM_Coverage_category,fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=I(4),alpha=I(0.8))+
xlab("LD Axis 1")+ylab("LD Axis 2")+theme_simple() +
theme(legend.title=element_blank())
#test significance of axis
anova(lm(FloristicSurveyLDAval$LD1~FloristicSurveysmall$GM_Coverage_category))anova(lm(FloristicSurveyLDAval$LD2~FloristicSurveysmall$GM_Coverage_category))Cluster analysis
# Ward Hierarchical Clustering with Bootstrapped p values
#do they cluster by population?
FloristicSurvey[is.na(FloristicSurvey)] <- 0
mydata<-t(data.frame(FloristicSurvey[,1],FloristicSurvey[,4:55]))
my.names <- mydata[1,]
colnames(mydata) <- my.names
mydata <- mydata[-1,]
d <- dist(t(mydata), method = "euclidean") # distance matrix
fit <- hclust(d, method = "ward.D")
fitd <- as.dendrogram(fit)
par(mfrow=c(3,2))
plot(fit,cex=0.4,main="Clustering with populations")
plot(cut(fitd, h=10)$upper, main="Upper tree of cut at h=10")
plot(cut(fitd, h=10)$lower[[1]], main="First branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[2]], main="Second branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[3]], main="Third branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[4]], main="Fourth branch of lower tree with cut at h=10")
#do they cluster by GM coverage?
FloristicSurvey[is.na(FloristicSurvey)] <- 0
#this will just make the dendogram easier to see
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("None", "ML", "MH", "High"))
mydata<-t(data.frame(FloristicSurvey[,56],FloristicSurvey[,4:55]))
my.names <- mydata[1,]
colnames(mydata) <- my.names
mydata <- mydata[-1,]
d <- dist(t(mydata), method = "euclidean") # distance matrix
fit <- hclust(d, method = "ward.D")
fitd <- as.dendrogram(fit)
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,5,30,50,70,95,100), labels = c("None", "Low", "Med_Low", "Med_High", "High", "Overtaken"))
par(mfrow=c(3,2))
plot(fit,cex=0.4,main="Clustering with populations")
plot(cut(fitd, h=10)$upper, main="Upper tree of cut at h=10")
plot(cut(fitd, h=10)$lower[[1]], main="First branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[2]], main="Second branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[3]], main="Third branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[4]], main="Fourth branch of lower tree with cut at h=10")
NMDS
# https://jonlefcheck.net/2012/10/24/nmds-tutorial-in-r/
# https://oliviarata.wordpress.com/2014/04/17/ordinations-in-ggplot2/
#Make a matrix with no row or column equal to 0 (do not enclude the env variable (GM COVERAGE))
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("Low", "Med_Low", "Med_High", "High"))
table(FloristicSurvey$GM_Coverage_category)##
## Low Med_Low Med_High High
## 82 24 58 34M <- as.matrix(FloristicSurvey[1:198,4:56])
M[is.na(M)] <- 0
rownames(M) <- FloristicSurvey$Population
M<-M[-c(130,133),]
M<-M[,-45]
GM_coverage_df <- data.frame(M[,52])
M<-M[,-52]
class(M) <- "numeric"
# with vegdist from Bray to seroeson: add binary = T
dist_FloristicSurvey <- vegdist(M, method = "bray", binary = T)
#The metaMDS analysis could have done the distance matrix internally but i would rather control it since i have presence/abscence
meta.nmds.FloristicSurvey2D <- metaMDS(dist_FloristicSurvey, k=2, trymax = 1000)## Run 0 stress 0.2097898
## Run 1 stress 0.2103351
## Run 2 stress 0.2099173
## ... Procrustes: rmse 0.01192537 max resid 0.1277007
## Run 3 stress 0.2126692
## Run 4 stress 0.2131916
## Run 5 stress 0.2130193
## Run 6 stress 0.2186913
## Run 7 stress 0.2098392
## ... Procrustes: rmse 0.009869736 max resid 0.1035989
## Run 8 stress 0.2352097
## Run 9 stress 0.211495
## Run 10 stress 0.2097224
## ... New best solution
## ... Procrustes: rmse 0.006331696 max resid 0.0529704
## Run 11 stress 0.2103191
## Run 12 stress 0.2122882
## Run 13 stress 0.2101326
## ... Procrustes: rmse 0.009866801 max resid 0.09557844
## Run 14 stress 0.2099405
## ... Procrustes: rmse 0.01634622 max resid 0.1525454
## Run 15 stress 0.2129822
## Run 16 stress 0.2102758
## Run 17 stress 0.2133365
## Run 18 stress 0.2102033
## ... Procrustes: rmse 0.006775596 max resid 0.07193472
## Run 19 stress 0.2097459
## ... Procrustes: rmse 0.008284148 max resid 0.09607869
## Run 20 stress 0.2117128
## Run 21 stress 0.2131278
## Run 22 stress 0.2099342
## ... Procrustes: rmse 0.01542433 max resid 0.152828
## Run 23 stress 0.2103332
## Run 24 stress 0.2102011
## ... Procrustes: rmse 0.006702283 max resid 0.07031149
## Run 25 stress 0.2202863
## Run 26 stress 0.2099359
## ... Procrustes: rmse 0.01602178 max resid 0.152552
## Run 27 stress 0.2098393
## ... Procrustes: rmse 0.01134596 max resid 0.1269412
## Run 28 stress 0.2146908
## Run 29 stress 0.2169865
## Run 30 stress 0.2129698
## Run 31 stress 0.2135891
## Run 32 stress 0.2170369
## Run 33 stress 0.2098652
## ... Procrustes: rmse 0.009361348 max resid 0.1212622
## Run 34 stress 0.2102694
## Run 35 stress 0.2159039
## Run 36 stress 0.2102638
## Run 37 stress 0.2163209
## Run 38 stress 0.2103664
## Run 39 stress 0.2126085
## Run 40 stress 0.2124092
## Run 41 stress 0.2129505
## Run 42 stress 0.2198644
## Run 43 stress 0.2133866
## Run 44 stress 0.2254004
## Run 45 stress 0.2192341
## Run 46 stress 0.213389
## Run 47 stress 0.2228373
## Run 48 stress 0.2099019
## ... Procrustes: rmse 0.01294875 max resid 0.1525675
## Run 49 stress 0.209759
## ... Procrustes: rmse 0.00652526 max resid 0.04372295
## Run 50 stress 0.2102149
## ... Procrustes: rmse 0.01700174 max resid 0.1532311
## Run 51 stress 0.228718
## Run 52 stress 0.2099745
## ... Procrustes: rmse 0.01454891 max resid 0.1533141
## Run 53 stress 0.2130074
## Run 54 stress 0.2126532
## Run 55 stress 0.2132173
## Run 56 stress 0.21023
## Run 57 stress 0.2136691
## Run 58 stress 0.2176221
## Run 59 stress 0.2125307
## Run 60 stress 0.2130349
## Run 61 stress 0.2101029
## ... Procrustes: rmse 0.01033814 max resid 0.09608037
## Run 62 stress 0.209758
## ... Procrustes: rmse 0.006722914 max resid 0.04813515
## Run 63 stress 0.2098131
## ... Procrustes: rmse 0.008581975 max resid 0.08773394
## Run 64 stress 0.2278306
## Run 65 stress 0.2102423
## Run 66 stress 0.2101894
## ... Procrustes: rmse 0.00970325 max resid 0.1026161
## Run 67 stress 0.212608
## Run 68 stress 0.2168897
## Run 69 stress 0.2099404
## ... Procrustes: rmse 0.0163191 max resid 0.1526328
## Run 70 stress 0.2131248
## Run 71 stress 0.2122657
## Run 72 stress 0.2133126
## Run 73 stress 0.2146229
## Run 74 stress 0.2099431
## ... Procrustes: rmse 0.01432381 max resid 0.1532026
## Run 75 stress 0.2173323
## Run 76 stress 0.2152484
## Run 77 stress 0.2117179
## Run 78 stress 0.2102342
## Run 79 stress 0.2098498
## ... Procrustes: rmse 0.008711821 max resid 0.0869956
## Run 80 stress 0.209855
## ... Procrustes: rmse 0.0111704 max resid 0.1073192
## Run 81 stress 0.2130124
## Run 82 stress 0.2099306
## ... Procrustes: rmse 0.01604905 max resid 0.1525495
## Run 83 stress 0.2148915
## Run 84 stress 0.2133183
## Run 85 stress 0.2152605
## Run 86 stress 0.2124853
## Run 87 stress 0.2112723
## Run 88 stress 0.2101393
## ... Procrustes: rmse 0.007157299 max resid 0.07225963
## Run 89 stress 0.212616
## Run 90 stress 0.2116955
## Run 91 stress 0.224748
## Run 92 stress 0.2097611
## ... Procrustes: rmse 0.004834871 max resid 0.02894467
## Run 93 stress 0.2134723
## Run 94 stress 0.2102688
## Run 95 stress 0.2098322
## ... Procrustes: rmse 0.0129102 max resid 0.1366845
## Run 96 stress 0.2099257
## ... Procrustes: rmse 0.01588553 max resid 0.1526738
## Run 97 stress 0.2131946
## Run 98 stress 0.2160033
## Run 99 stress 0.2103055
## Run 100 stress 0.2098845
## ... Procrustes: rmse 0.0126558 max resid 0.1525857
## Run 101 stress 0.2098822
## ... Procrustes: rmse 0.01203187 max resid 0.1530025
## Run 102 stress 0.209722
## ... New best solution
## ... Procrustes: rmse 0.005818339 max resid 0.0445689
## Run 103 stress 0.2133288
## Run 104 stress 0.2163187
## Run 105 stress 0.2161437
## Run 106 stress 0.2123118
## Run 107 stress 0.2129505
## Run 108 stress 0.2099305
## ... Procrustes: rmse 0.01562328 max resid 0.1533819
## Run 109 stress 0.2097666
## ... Procrustes: rmse 0.00932821 max resid 0.1214356
## Run 110 stress 0.2161913
## Run 111 stress 0.2101285
## ... Procrustes: rmse 0.007023751 max resid 0.07412396
## Run 112 stress 0.2097272
## ... Procrustes: rmse 0.00305532 max resid 0.02186644
## Run 113 stress 0.2118821
## Run 114 stress 0.2117277
## Run 115 stress 0.2098152
## ... Procrustes: rmse 0.01051393 max resid 0.09773606
## Run 116 stress 0.2133935
## Run 117 stress 0.2119973
## Run 118 stress 0.2129695
## Run 119 stress 0.2120935
## Run 120 stress 0.2134367
## Run 121 stress 0.2126096
## Run 122 stress 0.2129986
## Run 123 stress 0.2126035
## Run 124 stress 0.2097484
## ... Procrustes: rmse 0.002385626 max resid 0.02715788
## Run 125 stress 0.2370965
## Run 126 stress 0.2146537
## Run 127 stress 0.2099287
## ... Procrustes: rmse 0.01564831 max resid 0.1533717
## Run 128 stress 0.2098123
## ... Procrustes: rmse 0.01282204 max resid 0.1264842
## Run 129 stress 0.2150643
## Run 130 stress 0.2099154
## ... Procrustes: rmse 0.01450937 max resid 0.1534787
## Run 131 stress 0.2116439
## Run 132 stress 0.2100327
## ... Procrustes: rmse 0.01605158 max resid 0.1537851
## Run 133 stress 0.2115102
## Run 134 stress 0.2182735
## Run 135 stress 0.2103294
## Run 136 stress 0.2098103
## ... Procrustes: rmse 0.01094055 max resid 0.12632
## Run 137 stress 0.2166619
## Run 138 stress 0.2137642
## Run 139 stress 0.2126141
## Run 140 stress 0.2127359
## Run 141 stress 0.210165
## ... Procrustes: rmse 0.007980591 max resid 0.0748096
## Run 142 stress 0.2134329
## Run 143 stress 0.2159204
## Run 144 stress 0.2101602
## ... Procrustes: rmse 0.007766581 max resid 0.0746195
## Run 145 stress 0.2099195
## ... Procrustes: rmse 0.01187516 max resid 0.153661
## Run 146 stress 0.2154336
## Run 147 stress 0.2116886
## Run 148 stress 0.2101283
## ... Procrustes: rmse 0.007278672 max resid 0.07457893
## Run 149 stress 0.2097668
## ... Procrustes: rmse 0.003413671 max resid 0.03556374
## Run 150 stress 0.2097016
## ... New best solution
## ... Procrustes: rmse 0.004209457 max resid 0.02578612
## Run 151 stress 0.2101298
## ... Procrustes: rmse 0.0055622 max resid 0.07442458
## Run 152 stress 0.2112718
## Run 153 stress 0.2112806
## Run 154 stress 0.212974
## Run 155 stress 0.2101334
## ... Procrustes: rmse 0.007815866 max resid 0.07174455
## Run 156 stress 0.2097497
## ... Procrustes: rmse 0.004657534 max resid 0.02896196
## Run 157 stress 0.2137129
## Run 158 stress 0.2126415
## Run 159 stress 0.2097895
## ... Procrustes: rmse 0.005618096 max resid 0.02922783
## Run 160 stress 0.2097933
## ... Procrustes: rmse 0.009126156 max resid 0.1103368
## Run 161 stress 0.2097696
## ... Procrustes: rmse 0.006437338 max resid 0.04696808
## Run 162 stress 0.213439
## Run 163 stress 0.2103421
## Run 164 stress 0.2126159
## Run 165 stress 0.2102967
## Run 166 stress 0.2121353
## Run 167 stress 0.2224754
## Run 168 stress 0.2163261
## Run 169 stress 0.2248047
## Run 170 stress 0.210178
## ... Procrustes: rmse 0.01631992 max resid 0.153225
## Run 171 stress 0.2097509
## ... Procrustes: rmse 0.004694005 max resid 0.02844151
## Run 172 stress 0.2174803
## Run 173 stress 0.2101205
## ... Procrustes: rmse 0.005683145 max resid 0.07195479
## Run 174 stress 0.2150894
## Run 175 stress 0.2102499
## Run 176 stress 0.210227
## Run 177 stress 0.212503
## Run 178 stress 0.2147534
## Run 179 stress 0.2326758
## Run 180 stress 0.2097112
## ... Procrustes: rmse 0.005422487 max resid 0.04271508
## Run 181 stress 0.2116783
## Run 182 stress 0.2098008
## ... Procrustes: rmse 0.009201616 max resid 0.1118617
## Run 183 stress 0.2126199
## Run 184 stress 0.2099307
## ... Procrustes: rmse 0.01533 max resid 0.1530665
## Run 185 stress 0.2132181
## Run 186 stress 0.2098015
## ... Procrustes: rmse 0.009663798 max resid 0.1227011
## Run 187 stress 0.2112897
## Run 188 stress 0.2101601
## ... Procrustes: rmse 0.01772144 max resid 0.1664074
## Run 189 stress 0.210365
## Run 190 stress 0.2097791
## ... Procrustes: rmse 0.005172674 max resid 0.02794934
## Run 191 stress 0.2101423
## ... Procrustes: rmse 0.005995007 max resid 0.07402147
## Run 192 stress 0.2102844
## Run 193 stress 0.2098142
## ... Procrustes: rmse 0.009969578 max resid 0.126371
## Run 194 stress 0.2130184
## Run 195 stress 0.2103128
## Run 196 stress 0.2184558
## Run 197 stress 0.2102422
## Run 198 stress 0.213911
## Run 199 stress 0.2245082
## Run 200 stress 0.2150762
## Run 201 stress 0.2225779
## Run 202 stress 0.2129392
## Run 203 stress 0.2097176
## ... Procrustes: rmse 0.005990087 max resid 0.04087518
## Run 204 stress 0.2132634
## Run 205 stress 0.2098963
## ... Procrustes: rmse 0.01236197 max resid 0.1530079
## Run 206 stress 0.2126056
## Run 207 stress 0.224634
## Run 208 stress 0.2147059
## Run 209 stress 0.2198905
## Run 210 stress 0.213062
## Run 211 stress 0.209894
## ... Procrustes: rmse 0.01226484 max resid 0.1530733
## Run 212 stress 0.2247818
## Run 213 stress 0.209932
## ... Procrustes: rmse 0.01553434 max resid 0.1529955
## Run 214 stress 0.2136695
## Run 215 stress 0.2097927
## ... Procrustes: rmse 0.009248614 max resid 0.1160047
## Run 216 stress 0.2130952
## Run 217 stress 0.217344
## Run 218 stress 0.2125043
## Run 219 stress 0.2162334
## Run 220 stress 0.2098528
## ... Procrustes: rmse 0.01088873 max resid 0.09649557
## Run 221 stress 0.2099269
## ... Procrustes: rmse 0.01482054 max resid 0.1528368
## Run 222 stress 0.2144193
## Run 223 stress 0.2101275
## ... Procrustes: rmse 0.005750061 max resid 0.07361562
## Run 224 stress 0.212064
## Run 225 stress 0.2146626
## Run 226 stress 0.2137951
## Run 227 stress 0.2099296
## ... Procrustes: rmse 0.013993 max resid 0.1531029
## Run 228 stress 0.23434
## Run 229 stress 0.2154847
## Run 230 stress 0.2228331
## Run 231 stress 0.2144516
## Run 232 stress 0.2098424
## ... Procrustes: rmse 0.01039714 max resid 0.11132
## Run 233 stress 0.2097181
## ... Procrustes: rmse 0.00423394 max resid 0.02688357
## Run 234 stress 0.2145272
## Run 235 stress 0.2144672
## Run 236 stress 0.2120379
## Run 237 stress 0.2161529
## Run 238 stress 0.2117142
## Run 239 stress 0.2102623
## Run 240 stress 0.2161653
## Run 241 stress 0.2099349
## ... Procrustes: rmse 0.01458496 max resid 0.1531812
## Run 242 stress 0.2125044
## Run 243 stress 0.2126944
## Run 244 stress 0.2258701
## Run 245 stress 0.2267006
## Run 246 stress 0.2134065
## Run 247 stress 0.2129694
## Run 248 stress 0.2145001
## Run 249 stress 0.2123909
## Run 250 stress 0.2130058
## Run 251 stress 0.2113364
## Run 252 stress 0.2099315
## ... Procrustes: rmse 0.0153496 max resid 0.1530516
## Run 253 stress 0.2163246
## Run 254 stress 0.2116937
## Run 255 stress 0.2099188
## ... Procrustes: rmse 0.0150663 max resid 0.1530329
## Run 256 stress 0.2103265
## Run 257 stress 0.2102338
## Run 258 stress 0.2097157
## ... Procrustes: rmse 0.008446499 max resid 0.0961253
## Run 259 stress 0.2131688
## Run 260 stress 0.2099365
## ... Procrustes: rmse 0.01558455 max resid 0.1529922
## Run 261 stress 0.2101178
## ... Procrustes: rmse 0.005679899 max resid 0.07294909
## Run 262 stress 0.2126052
## Run 263 stress 0.2101146
## ... Procrustes: rmse 0.005736433 max resid 0.07349935
## Run 264 stress 0.2097854
## ... Procrustes: rmse 0.009056932 max resid 0.1018646
## Run 265 stress 0.2097486
## ... Procrustes: rmse 0.004576749 max resid 0.02854606
## Run 266 stress 0.2181454
## Run 267 stress 0.2101606
## ... Procrustes: rmse 0.01047877 max resid 0.1044486
## Run 268 stress 0.2099258
## ... Procrustes: rmse 0.01433792 max resid 0.1530823
## Run 269 stress 0.211862
## Run 270 stress 0.2097949
## ... Procrustes: rmse 0.006538335 max resid 0.04690468
## Run 271 stress 0.2102378
## Run 272 stress 0.2148592
## Run 273 stress 0.212043
## Run 274 stress 0.2148504
## Run 275 stress 0.2116562
## Run 276 stress 0.2101112
## ... Procrustes: rmse 0.01043439 max resid 0.1013377
## Run 277 stress 0.2130211
## Run 278 stress 0.2103682
## Run 279 stress 0.2099322
## ... Procrustes: rmse 0.0154728 max resid 0.1529921
## Run 280 stress 0.2098838
## ... Procrustes: rmse 0.01184859 max resid 0.1531203
## Run 281 stress 0.2148362
## Run 282 stress 0.2102986
## Run 283 stress 0.2126003
## Run 284 stress 0.2129711
## Run 285 stress 0.2098408
## ... Procrustes: rmse 0.00992117 max resid 0.122135
## Run 286 stress 0.2131976
## Run 287 stress 0.2188947
## Run 288 stress 0.2125095
## Run 289 stress 0.2125251
## Run 290 stress 0.2213913
## Run 291 stress 0.2103218
## Run 292 stress 0.2103556
## Run 293 stress 0.2102009
## ... Procrustes: rmse 0.006264532 max resid 0.07087182
## Run 294 stress 0.2128862
## Run 295 stress 0.2293045
## Run 296 stress 0.2101116
## ... Procrustes: rmse 0.005800042 max resid 0.0737227
## Run 297 stress 0.210243
## Run 298 stress 0.2134944
## Run 299 stress 0.2172719
## Run 300 stress 0.2162612
## Run 301 stress 0.2434299
## Run 302 stress 0.2102694
## Run 303 stress 0.2134215
## Run 304 stress 0.2127281
## Run 305 stress 0.2099516
## ... Procrustes: rmse 0.01221039 max resid 0.1536921
## Run 306 stress 0.2178865
## Run 307 stress 0.2097583
## ... Procrustes: rmse 0.005003836 max resid 0.02866227
## Run 308 stress 0.2117147
## Run 309 stress 0.2098077
## ... Procrustes: rmse 0.009842923 max resid 0.1219967
## Run 310 stress 0.2165044
## Run 311 stress 0.212987
## Run 312 stress 0.2103349
## Run 313 stress 0.2099146
## ... Procrustes: rmse 0.01485102 max resid 0.1530597
## Run 314 stress 0.2169177
## Run 315 stress 0.2247716
## Run 316 stress 0.2098862
## ... Procrustes: rmse 0.01220158 max resid 0.1531984
## Run 317 stress 0.2133246
## Run 318 stress 0.2167595
## Run 319 stress 0.2098052
## ... Procrustes: rmse 0.00953829 max resid 0.1244171
## Run 320 stress 0.2112655
## Run 321 stress 0.2116529
## Run 322 stress 0.2225972
## Run 323 stress 0.2102479
## Run 324 stress 0.2113148
## Run 325 stress 0.2175125
## Run 326 stress 0.2120987
## Run 327 stress 0.210256
## Run 328 stress 0.2134304
## Run 329 stress 0.2102517
## Run 330 stress 0.2129867
## Run 331 stress 0.2150343
## Run 332 stress 0.2098352
## ... Procrustes: rmse 0.00919755 max resid 0.09367652
## Run 333 stress 0.2168579
## Run 334 stress 0.2118603
## Run 335 stress 0.2101493
## ... Procrustes: rmse 0.01068699 max resid 0.1049335
## Run 336 stress 0.2102338
## Run 337 stress 0.2099255
## ... Procrustes: rmse 0.01505611 max resid 0.1528571
## Run 338 stress 0.2099145
## ... Procrustes: rmse 0.0116063 max resid 0.0973198
## Run 339 stress 0.2099309
## ... Procrustes: rmse 0.01533235 max resid 0.1530627
## Run 340 stress 0.2097275
## ... Procrustes: rmse 0.00395705 max resid 0.02800491
## Run 341 stress 0.2153898
## Run 342 stress 0.2394691
## Run 343 stress 0.21273
## Run 344 stress 0.2102831
## Run 345 stress 0.2097241
## ... Procrustes: rmse 0.004416164 max resid 0.0280246
## Run 346 stress 0.2098928
## ... Procrustes: rmse 0.01240304 max resid 0.152975
## Run 347 stress 0.2097859
## ... Procrustes: rmse 0.005706258 max resid 0.06201415
## Run 348 stress 0.2134563
## Run 349 stress 0.2121168
## Run 350 stress 0.210378
## Run 351 stress 0.2103764
## Run 352 stress 0.2130106
## Run 353 stress 0.2170193
## Run 354 stress 0.2098083
## ... Procrustes: rmse 0.005275467 max resid 0.02941205
## Run 355 stress 0.2119085
## Run 356 stress 0.2103907
## Run 357 stress 0.2102628
## Run 358 stress 0.2101609
## ... Procrustes: rmse 0.005813887 max resid 0.07415652
## Run 359 stress 0.2099309
## ... Procrustes: rmse 0.01456246 max resid 0.1532625
## Run 360 stress 0.2175221
## Run 361 stress 0.2217646
## Run 362 stress 0.2102498
## Run 363 stress 0.2102884
## Run 364 stress 0.2169011
## Run 365 stress 0.2132479
## Run 366 stress 0.2099212
## ... Procrustes: rmse 0.01513998 max resid 0.1529961
## Run 367 stress 0.2131867
## Run 368 stress 0.2131265
## Run 369 stress 0.2116568
## Run 370 stress 0.2231066
## Run 371 stress 0.2098369
## ... Procrustes: rmse 0.009911654 max resid 0.1142778
## Run 372 stress 0.2097854
## ... Procrustes: rmse 0.008789799 max resid 0.1177037
## Run 373 stress 0.2097412
## ... Procrustes: rmse 0.007941668 max resid 0.08885918
## Run 374 stress 0.2102052
## Run 375 stress 0.2097891
## ... Procrustes: rmse 0.006550043 max resid 0.05198666
## Run 376 stress 0.213326
## Run 377 stress 0.2097025
## ... Procrustes: rmse 0.002876355 max resid 0.01595654
## Run 378 stress 0.2099077
## ... Procrustes: rmse 0.01204289 max resid 0.1534777
## Run 379 stress 0.212503
## Run 380 stress 0.2148804
## Run 381 stress 0.2131164
## Run 382 stress 0.2141244
## Run 383 stress 0.2098294
## ... Procrustes: rmse 0.01036347 max resid 0.1322424
## Run 384 stress 0.2102495
## Run 385 stress 0.2102679
## Run 386 stress 0.2179382
## Run 387 stress 0.2198484
## Run 388 stress 0.2097184
## ... Procrustes: rmse 0.004210321 max resid 0.02849628
## Run 389 stress 0.212339
## Run 390 stress 0.2148237
## Run 391 stress 0.2099264
## ... Procrustes: rmse 0.01534003 max resid 0.1529914
## Run 392 stress 0.2102592
## Run 393 stress 0.2101307
## ... Procrustes: rmse 0.005484187 max resid 0.07336651
## Run 394 stress 0.209805
## ... Procrustes: rmse 0.009558297 max resid 0.120224
## Run 395 stress 0.2163558
## Run 396 stress 0.2102291
## Run 397 stress 0.2161451
## Run 398 stress 0.2126068
## Run 399 stress 0.2099261
## ... Procrustes: rmse 0.01509752 max resid 0.1528582
## Run 400 stress 0.2347234
## Run 401 stress 0.2112807
## Run 402 stress 0.2097053
## ... Procrustes: rmse 0.004117119 max resid 0.02287897
## Run 403 stress 0.2107541
## Run 404 stress 0.209775
## ... Procrustes: rmse 0.009373035 max resid 0.1202368
## Run 405 stress 0.2151085
## Run 406 stress 0.237914
## Run 407 stress 0.2131815
## Run 408 stress 0.213702
## Run 409 stress 0.2277012
## Run 410 stress 0.2097137
## ... Procrustes: rmse 0.003256294 max resid 0.02391741
## Run 411 stress 0.2126363
## Run 412 stress 0.2098564
## ... Procrustes: rmse 0.009359887 max resid 0.09464593
## Run 413 stress 0.2096806
## ... New best solution
## ... Procrustes: rmse 0.003771369 max resid 0.02563352
## Run 414 stress 0.2101414
## ... Procrustes: rmse 0.01169878 max resid 0.1053441
## Run 415 stress 0.211258
## Run 416 stress 0.21295
## Run 417 stress 0.2146358
## Run 418 stress 0.2099262
## ... Procrustes: rmse 0.01465429 max resid 0.153668
## Run 419 stress 0.2097208
## ... Procrustes: rmse 0.009068168 max resid 0.1005321
## Run 420 stress 0.2099474
## ... Procrustes: rmse 0.01508528 max resid 0.1535732
## Run 421 stress 0.2097157
## ... Procrustes: rmse 0.004738058 max resid 0.04168663
## Run 422 stress 0.2117891
## Run 423 stress 0.2179079
## Run 424 stress 0.2099126
## ... Procrustes: rmse 0.01355478 max resid 0.1535678
## Run 425 stress 0.2102493
## Run 426 stress 0.2220886
## Run 427 stress 0.2099387
## ... Procrustes: rmse 0.01348499 max resid 0.1537097
## Run 428 stress 0.2098314
## ... Procrustes: rmse 0.0076032 max resid 0.09563233
## Run 429 stress 0.2112498
## Run 430 stress 0.2101959
## Run 431 stress 0.2134495
## Run 432 stress 0.2134059
## Run 433 stress 0.2098668
## ... Procrustes: rmse 0.0103811 max resid 0.08894867
## Run 434 stress 0.2226698
## Run 435 stress 0.2130046
## Run 436 stress 0.2098851
## ... Procrustes: rmse 0.01242216 max resid 0.1534893
## Run 437 stress 0.2157712
## Run 438 stress 0.2166472
## Run 439 stress 0.2180417
## Run 440 stress 0.2097523
## ... Procrustes: rmse 0.003115454 max resid 0.02803743
## Run 441 stress 0.2148538
## Run 442 stress 0.2103
## Run 443 stress 0.2117119
## Run 444 stress 0.2097551
## ... Procrustes: rmse 0.006921319 max resid 0.08624046
## Run 445 stress 0.2099123
## ... Procrustes: rmse 0.0138407 max resid 0.1535531
## Run 446 stress 0.2097192
## ... Procrustes: rmse 0.005393117 max resid 0.04242808
## Run 447 stress 0.2148501
## Run 448 stress 0.2195586
## Run 449 stress 0.2096752
## ... New best solution
## ... Procrustes: rmse 0.00120412 max resid 0.01553874
## Run 450 stress 0.2099407
## ... Procrustes: rmse 0.01512919 max resid 0.1536836
## Run 451 stress 0.2177598
## Run 452 stress 0.2133447
## Run 453 stress 0.2097955
## ... Procrustes: rmse 0.003885868 max resid 0.02744706
## Run 454 stress 0.2097694
## ... Procrustes: rmse 0.005719434 max resid 0.04377941
## Run 455 stress 0.2101171
## ... Procrustes: rmse 0.01020592 max resid 0.09090328
## Run 456 stress 0.209809
## ... Procrustes: rmse 0.004514653 max resid 0.02835068
## Run 457 stress 0.2102806
## Run 458 stress 0.2099179
## ... Procrustes: rmse 0.01418331 max resid 0.1533645
## Run 459 stress 0.2126038
## Run 460 stress 0.2101389
## ... Procrustes: rmse 0.008306194 max resid 0.07348261
## Run 461 stress 0.2220578
## Run 462 stress 0.2096898
## ... Procrustes: rmse 0.002116752 max resid 0.010882
## Run 463 stress 0.2120487
## Run 464 stress 0.2103361
## Run 465 stress 0.2194496
## Run 466 stress 0.2175507
## Run 467 stress 0.2098872
## ... Procrustes: rmse 0.0119993 max resid 0.1534727
## Run 468 stress 0.2132519
## Run 469 stress 0.2129907
## Run 470 stress 0.2102533
## Run 471 stress 0.212572
## Run 472 stress 0.2102243
## Run 473 stress 0.2112709
## Run 474 stress 0.2349923
## Run 475 stress 0.2126826
## Run 476 stress 0.2154687
## Run 477 stress 0.2112593
## Run 478 stress 0.2102025
## Run 479 stress 0.2147848
## Run 480 stress 0.2097926
## ... Procrustes: rmse 0.005780708 max resid 0.04216166
## Run 481 stress 0.2146509
## Run 482 stress 0.2260598
## Run 483 stress 0.2116651
## Run 484 stress 0.2098007
## ... Procrustes: rmse 0.004470413 max resid 0.02955439
## Run 485 stress 0.2099029
## ... Procrustes: rmse 0.01242719 max resid 0.1534549
## Run 486 stress 0.2102856
## Run 487 stress 0.2116572
## Run 488 stress 0.2124048
## Run 489 stress 0.211685
## Run 490 stress 0.2096751
## ... New best solution
## ... Procrustes: rmse 0.0005534366 max resid 0.004358859
## ... Similar to previous best
## *** Solution reachedstr(meta.nmds.FloristicSurvey2D)## List of 35
## $ nobj : int 196
## $ nfix : int 0
## $ ndim : num 2
## $ ndis : int 19110
## $ ngrp : int 1
## $ diss : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
## $ iidx : int [1:19110] 117 111 67 129 186 123 126 129 129 170 ...
## $ jidx : int [1:19110] 106 103 2 122 128 121 122 127 126 3 ...
## $ xinit : num [1:392] 0.849 0.803 0.633 0.357 0.187 ...
## $ istart : int 1
## $ isform : int 1
## $ ities : int 1
## $ iregn : int 1
## $ iscal : int 1
## $ maxits : int 200
## $ sratmx : num 1
## $ strmin : num 1e-04
## $ sfgrmn : num 1e-07
## $ dist : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
## $ dhat : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
## $ points : num [1:196, 1:2] -0.136 -0.149 -0.154 -0.105 -0.103 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:196] "7" "7" "7" "7" ...
## .. ..$ : chr [1:2] "MDS1" "MDS2"
## ..- attr(*, "centre")= logi TRUE
## ..- attr(*, "pc")= logi TRUE
## ..- attr(*, "halfchange")= logi FALSE
## ..- attr(*, "internalscaling")= num 4.31
## $ stress : num 0.21
## $ grstress : num 0.21
## $ iters : int 148
## $ icause : int 3
## $ call : language metaMDS(comm = dist_FloristicSurvey, k = 2, trymax = 1000)
## $ model : chr "global"
## $ distmethod: chr "binary bray"
## $ distcall : chr "vegdist(x = M, method = \"bray\", binary = T)"
## $ distance : chr "binary bray"
## $ converged : logi TRUE
## $ tries : num 490
## $ engine : chr "monoMDS"
## $ species : logi NA
## $ data : chr "dist_FloristicSurvey"
## - attr(*, "class")= chr [1:2] "metaMDS" "monoMDS"stressplot(meta.nmds.FloristicSurvey2D)
FloristicSurvey_envfit <- envfit(meta.nmds.FloristicSurvey2D, env = GM_coverage_df, perm = 999) #standard envfit
FloristicSurvey_envfit##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## M...52.High 0.0453 0.0401
## M...52.Low -0.0896 -0.0251
## M...52.Med_High 0.0782 -0.0018
## M...52.Med_Low 0.0455 0.0313
##
## Goodness of fit:
## r2 Pr(>r)
## M...52. 0.1181 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999#data for plotting
##NMDS points
FloristicSurvey.NMDS.data<-GM_coverage_df
FloristicSurvey.NMDS.data$NMDS1<-meta.nmds.FloristicSurvey2D$points[ ,1]
FloristicSurvey.NMDS.data$NMDS2<-meta.nmds.FloristicSurvey2D$points[ ,2]
colnames(FloristicSurvey.NMDS.data)[1] <- "GM_Coverage"
# data for the envfit arrows
env.scores.FloristicSurvey <- as.data.frame(scores(FloristicSurvey_envfit, display = "vectors")) #extracts relevant scores from envifit
env.scores.FloristicSurvey <- cbind(env.scores.FloristicSurvey, env.variables = rownames(env.scores.FloristicSurvey)) #and then gives them their names
# function for ellipsess - just run this, is used later
#taken from the excellent stackoverflow Q+A: http://stackoverflow.com/questions/13794419/plotting-ordiellipse-function-from-vegan-package-onto-nmds-plot-created-in-ggplo
veganCovEllipse <- function (cov, center = c(0, 0), scale = 1, npoints = 100)
{
theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))
}
#data for ellipse, use GM coverage
df_ell.FSurvey.GM_coverage <- data.frame() #sets up a data frame before running the function.
for(g in levels(FloristicSurvey.NMDS.data$GM_Coverage)){
df_ell.FSurvey.GM_coverage <- rbind(df_ell.FSurvey.GM_coverage, cbind(as.data.frame(with(FloristicSurvey.NMDS.data [FloristicSurvey.NMDS.data$GM_Coverage==g,],
veganCovEllipse(cov.wt(cbind(NMDS1,NMDS2),wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS2)))))
,GM_coverage=g))
}
# data for labelling the ellipse
NMDS.mean.FloristicSurvey=aggregate(FloristicSurvey.NMDS.data[ ,c("NMDS1", "NMDS2")],
list(group = FloristicSurvey.NMDS.data$GM_Coverage), mean)
## finally plotting.
mult <- 1 #multiplier for the arrows and text for envfit below. You can change this and then rerun the plot command.
ggplot(data = FloristicSurvey.NMDS.data, aes(y = NMDS2, x = NMDS1))+
geom_path(data = df_ell.FSurvey.GM_coverage, aes(x = NMDS1, y = NMDS2, group = df_ell.FSurvey.GM_coverage$GM_coverage, color=df_ell.FSurvey.GM_coverage$GM_coverage))+
geom_point( aes(color = FloristicSurvey.NMDS.data$GM_Coverage), size = 1)+
scale_color_manual(values = c(14,5,19,2))+
coord_cartesian(xlim = c(-1,1.5))+
theme_simple()+ guides(color=guide_legend("Garlic Mustard Coverage"))
PCA
FloristicSurvey<-FloristicSurvey[,-56]
FloristicSurvey.pca <- prcomp(FloristicSurvey[,4:(ncol(FloristicSurvey))],center = TRUE)
summary(FloristicSurvey.pca)## Importance of components%s:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 0.6406 0.5639 0.49314 0.45065 0.42882 0.3954
## Proportion of Variance 0.1365 0.1058 0.08088 0.06755 0.06116 0.0520
## Cumulative Proportion 0.1365 0.2422 0.32311 0.39066 0.45182 0.5038
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.37368 0.3608 0.34564 0.33141 0.30059 0.27972
## Proportion of Variance 0.04644 0.0433 0.03974 0.03653 0.03005 0.02602
## Cumulative Proportion 0.55027 0.5936 0.63330 0.66983 0.69989 0.72591
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 0.26731 0.23707 0.23511 0.23170 0.23140 0.2159
## Proportion of Variance 0.02377 0.01869 0.01839 0.01786 0.01781 0.0155
## Cumulative Proportion 0.74968 0.76837 0.78675 0.80461 0.82242 0.8379
## PC19 PC20 PC21 PC22 PC23 PC24
## Standard deviation 0.20362 0.1992 0.1835 0.17727 0.17130 0.16138
## Proportion of Variance 0.01379 0.0132 0.0112 0.01045 0.00976 0.00866
## Cumulative Proportion 0.85171 0.8649 0.8761 0.88656 0.89632 0.90499
## PC25 PC26 PC27 PC28 PC29 PC30
## Standard deviation 0.16000 0.15304 0.14850 0.14015 0.13172 0.12966
## Proportion of Variance 0.00851 0.00779 0.00733 0.00653 0.00577 0.00559
## Cumulative Proportion 0.91350 0.92129 0.92863 0.93516 0.94093 0.94652
## PC31 PC32 PC33 PC34 PC35 PC36
## Standard deviation 0.12335 0.11774 0.11421 0.11329 0.10677 0.10455
## Proportion of Variance 0.00506 0.00461 0.00434 0.00427 0.00379 0.00364
## Cumulative Proportion 0.95158 0.95619 0.96053 0.96480 0.96859 0.97223
## PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 0.10160 0.09716 0.09407 0.08205 0.07828 0.07736
## Proportion of Variance 0.00343 0.00314 0.00294 0.00224 0.00204 0.00199
## Cumulative Proportion 0.97566 0.97880 0.98174 0.98398 0.98602 0.98801
## PC43 PC44 PC45 PC46 PC47 PC48
## Standard deviation 0.07252 0.06984 0.06822 0.06603 0.06268 0.05955
## Proportion of Variance 0.00175 0.00162 0.00155 0.00145 0.00131 0.00118
## Cumulative Proportion 0.98976 0.99138 0.99293 0.99438 0.99569 0.99687
## PC49 PC50 PC51 PC52
## Standard deviation 0.05928 0.05691 0.05170 4.972e-17
## Proportion of Variance 0.00117 0.00108 0.00089 0.000e+00
## Cumulative Proportion 0.99803 0.99911 1.00000 1.000e+00fviz_eig(FloristicSurvey.pca,addlabels=TRUE,hjust = -0.3)+theme_simple()
round(FloristicSurvey.pca$sdev^2/sum(FloristicSurvey.pca$sdev^2)*100,2)## [1] 13.65 10.58 8.09 6.75 6.12 5.20 4.64 4.33 3.97 3.65 3.01
## [12] 2.60 2.38 1.87 1.84 1.79 1.78 1.55 1.38 1.32 1.12 1.05
## [23] 0.98 0.87 0.85 0.78 0.73 0.65 0.58 0.56 0.51 0.46 0.43
## [34] 0.43 0.38 0.36 0.34 0.31 0.29 0.22 0.20 0.20 0.17 0.16
## [45] 0.15 0.15 0.13 0.12 0.12 0.11 0.09 0.00# About 25% explained by my first 2 eigenvectors. 50% by the first 6.Modeling (redundancy analysis)
species comp ~ GMcov * pop + soil characteristics inside and outside as covariates?
species<-names(FloristicSurvey[4:(length(FloristicSurvey))])
SpeciesGMCorrelation <- data.frame(matrix(ncol = 3, nrow = 0))
names(SpeciesGMCorrelation)<-c("Species", "PValue","signif")
for(i in 4:(length(FloristicSurvey))){
SpeciesGMCorrelation[i-3,1]<-names(FloristicSurvey[i])
SpeciesGMCorrelation[i-3,2]<-summary(glm(FloristicSurvey[,i]~GM_Coverage, data=FloristicSurvey, family = binomial))$coefficients[8]
}## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not convergefor(i in 1:length(species)) if( SpeciesGMCorrelation[i,2] <= 0.05){
SpeciesGMCorrelation[i,3]<-c("TRUE")
}else{
SpeciesGMCorrelation[i,3]<-c("FALSE")
}
length(SpeciesGMCorrelation$signif[SpeciesGMCorrelation[,3]==TRUE])/length(SpeciesGMCorrelation$signif)*100## [1] 25SpeciesGMCorrelation$Species[SpeciesGMCorrelation[,3]==TRUE]## [1] "Acer_Saccharum" "Galium_circaezans"
## [3] "Ulmus_americana" "Aliaria_petiolata"
## [5] "Parthenocissus_quinquefolia" "Rhubarb"
## [7] "Trifolium_sp" "Rubus_occidentalis"
## [9] "Circeae_lutetiana" "Clematis_virginiana"
## [11] "Rhamnus_cathartica" "Unknown_8"
## [13] "unknown_12"# Redundancy analysis
sp.comp <- FloristicSurvey[,4:(length(FloristicSurvey))]
var <- FloristicSurvey[,1:3]
var$Quadrat <- as.character(var$Quadrat)
for(i in 1:length(var$Quadrat))if(var$Population[i] == "7"){
var$Quadrat[i]<-paste("7",var$Quadrat[i],sep="")
}else if(var$Population[i] == "3"){
var$Quadrat[i]<-paste("3",var$Quadrat[i],sep="")
}else if(var$Population[i] == "14"){
var$Quadrat[i]<-paste("14",var$Quadrat[i],sep="")
}else if(var$Population[i] == "13"){
var$Quadrat[i]<-paste("13",var$Quadrat[i],sep="")
}else if(var$Population[i] == "1"){
var$Quadrat[i]<-paste("1",var$Quadrat[i],sep="")
}
var$Population <- as.factor(var$Population)
formulaRDA <- rda(sp.comp ~ GM_Coverage + Condition(Population) , data=var)
formulaRDA <- rda(sp.comp ~ GM_Coverage , data=var)
head(summary(formulaRDA))##
## Call:
## rda(formula = sp.comp ~ GM_Coverage, data = var)
##
## Partitioning of variance:
## Inertia Proportion
## Total 3.0066 1.00000
## Constrained 0.1451 0.04825
## Unconstrained 2.8615 0.95175
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 PC1 PC2 PC3 PC4 PC5
## Eigenvalue 0.14507 0.3830 0.25896 0.23497 0.19681 0.18378
## Proportion Explained 0.04825 0.1274 0.08613 0.07815 0.06546 0.06113
## Cumulative Proportion 0.04825 0.1757 0.26178 0.33993 0.40539 0.46652
## PC6 PC7 PC8 PC9 PC10 PC11
## Eigenvalue 0.14105 0.1302 0.12555 0.11261 0.10979 0.09025
## Proportion Explained 0.04691 0.0433 0.04176 0.03746 0.03652 0.03002
## Cumulative Proportion 0.51343 0.5567 0.59849 0.63595 0.67246 0.70248
## PC12 PC13 PC14 PC15 PC16 PC17
## Eigenvalue 0.07818 0.07134 0.05604 0.05392 0.05355 0.05152
## Proportion Explained 0.02600 0.02373 0.01864 0.01793 0.01781 0.01714
## Cumulative Proportion 0.72848 0.75221 0.77085 0.78878 0.80659 0.82373
## PC18 PC19 PC20 PC21 PC22 PC23
## Eigenvalue 0.04588 0.04146 0.0394 0.03355 0.03051 0.02929
## Proportion Explained 0.01526 0.01379 0.0131 0.01116 0.01015 0.00974
## Cumulative Proportion 0.83899 0.85278 0.8659 0.87704 0.88719 0.89693
## PC24 PC25 PC26 PC27 PC28 PC29
## Eigenvalue 0.02567 0.02559 0.02324 0.02183 0.01964 0.01731
## Proportion Explained 0.00854 0.00851 0.00773 0.00726 0.00653 0.00576
## Cumulative Proportion 0.90547 0.91398 0.92171 0.92897 0.93551 0.94126
## PC30 PC31 PC32 PC33 PC34 PC35
## Eigenvalue 0.01681 0.01508 0.01386 0.01304 0.01283 0.01139
## Proportion Explained 0.00559 0.00501 0.00461 0.00434 0.00427 0.00379
## Cumulative Proportion 0.94686 0.95187 0.95648 0.96082 0.96508 0.96887
## PC36 PC37 PC38 PC39 PC40 PC41
## Eigenvalue 0.01089 0.01027 0.009383 0.008708 0.006676 0.006114
## Proportion Explained 0.00362 0.00342 0.003120 0.002900 0.002220 0.002030
## Cumulative Proportion 0.97249 0.97591 0.979030 0.981930 0.984150 0.986180
## PC42 PC43 PC44 PC45 PC46
## Eigenvalue 0.005879 0.005254 0.004851 0.004653 0.004186
## Proportion Explained 0.001960 0.001750 0.001610 0.001550 0.001390
## Cumulative Proportion 0.988140 0.989880 0.991500 0.993050 0.994440
## PC47 PC48 PC49 PC50 PC51
## Eigenvalue 0.003919 0.003546 0.003458 0.003175 0.002623
## Proportion Explained 0.001300 0.001180 0.001150 0.001060 0.000870
## Cumulative Proportion 0.995740 0.996920 0.998070 0.999130 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1
## Eigenvalue 0.1451
## Proportion Explained 1.0000
## Cumulative Proportion 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 4.933262
##
##
## Species scores
##
## RDA1 PC1 PC2 PC3 PC4 PC5
## Galium_aparine -0.03884 0.95611 -0.68111 0.398026 0.56161 -0.16776
## Acer_Saccharum -0.39518 -0.59225 -0.76442 -0.724167 0.39649 0.40523
## Unknown_Grass -0.09710 0.83165 0.29277 -0.886732 -0.01378 -0.05976
## Galium_circaezans -0.09565 0.08758 -0.06161 -0.009061 0.03333 -0.06610
## Ulmus_americana -0.18939 0.19041 -0.59667 -0.292853 -0.51085 -0.72307
## False_Solomon_Seal -0.05895 0.05660 0.00144 -0.025574 0.01268 -0.01885
## ....
##
##
## Site scores (weighted sums of species scores)
##
## RDA1 PC1 PC2 PC3 PC4 PC5
## row1 -0.8792 0.3117 -0.4541 -0.3298 0.088931 -0.59497
## row2 -0.9951 0.2201 -0.4529 -0.3803 -0.087225 -0.39023
## row3 -0.6743 0.1376 0.2111 -0.5348 -0.312633 0.51839
## row4 -0.8687 0.5073 -0.3985 -0.4795 -0.280767 -0.08044
## row5 -0.7576 0.4642 -0.1132 -0.3390 0.047862 0.39295
## row6 -0.5475 0.3770 0.2156 0.1561 0.003758 -0.18769
## ....
##
##
## Site constraints (linear combinations of constraining variables)
##
## RDA1 PC1 PC2 PC3 PC4 PC5
## row1 -0.4053 0.3117 -0.4541 -0.3298 0.088931 -0.59497
## row2 -0.4053 0.2201 -0.4529 -0.3803 -0.087225 -0.39023
## row3 -0.4591 0.1376 0.2111 -0.5348 -0.312633 0.51839
## row4 -0.4591 0.5073 -0.3985 -0.4795 -0.280767 -0.08044
## row5 -0.4591 0.4642 -0.1132 -0.3390 0.047862 0.39295
## row6 -0.4591 0.3770 0.2156 0.1561 0.003758 -0.18769
## ....
##
##
## Biplot scores for constraining variables
##
## RDA1 PC1 PC2 PC3 PC4 PC5
## GM_Coverage 1 0 0 0 0 0RsquareAdj(formulaRDA)## $r.squared
## [1] 0.04825094
##
## $adj.r.squared
## [1] 0.04339508plot(formulaRDA, scaling = 2)
anova.cca(formulaRDA,step=1000)anova.cca(formulaRDA,by="axis", step=1000)*Note that I redid this with GM_coverage_catogory and i didnt find ANY significance!
Soil Characteristics
Soil_characteristics$Location<-as.factor(Soil_characteristics$Location)
GLmodel_soil <- glm(Location ~Agg_stability*C_percent*N_percent,family=binomial(link='logit'),data=Soil_characteristics)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(GLmodel_soil)##
## Call:
## glm(formula = Location ~ Agg_stability * C_percent * N_percent,
## family = binomial(link = "logit"), data = Soil_characteristics)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5019 -0.9103 -0.1493 0.7927 2.1777
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -40.68 40.25 -1.011 0.312
## Agg_stability 57.36 62.43 0.919 0.358
## C_percent 28.40 18.23 1.558 0.119
## N_percent -178.24 138.30 -1.289 0.197
## Agg_stability:C_percent -40.73 26.92 -1.513 0.130
## Agg_stability:N_percent 253.80 211.09 1.202 0.229
## C_percent:N_percent -14.38 10.04 -1.432 0.152
## Agg_stability:C_percent:N_percent 20.95 15.18 1.380 0.168
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 18.370 on 12 degrees of freedom
## AIC: 34.37
##
## Number of Fisher Scoring iterations: 8Soil_characteristics2 <- data.frame(Soil_characteristics$Population, Soil_characteristics$Location, log(Soil_characteristics$pH), log(Soil_characteristics$Agg_stability), log(Soil_characteristics$C_percent), log(Soil_characteristics$N_percent))
Soil_characteristics2<-rename(Soil_characteristics2, c("Soil_characteristics.Population"="Population", "Soil_characteristics.Location"="Location", "log.Soil_characteristics.pH."="pH", "log.Soil_characteristics.Agg_stability."="Agg_stability", "log.Soil_characteristics.C_percent."="Percent","log.Soil_characteristics.N_percent."="Percent" ))## Warning: The plyr::rename operation has created duplicates for the
## following name(s): (`Percent`)df2 <- reshape2::melt(Soil_characteristics2, id.var=c("Population","Location"))
df2 <- rename(df2, c("value"="Log_value"))
par(mfrow=c(2,2))
ggplot(Soil_characteristics, aes(Location,Agg_stability,fill=Location))+ geom_boxplot() +
theme_simple()
ggplot(Soil_characteristics, aes(Location,pH,fill=Location))+ geom_boxplot() +
theme_simple()
ggplot(Soil_characteristics, aes(Location,C_percent,fill=Location))+ geom_boxplot() +
theme_simple()
ggplot(Soil_characteristics, aes(Location,N_percent,fill=Location))+ geom_boxplot() +
theme_simple()
Soil_characteristics$Population<-as.factor(Soil_characteristics$Population)
GLMER_soil <- glmer(Location~scale(Agg_stability)*scale(C_percent)*scale(N_percent)+(1|Population),family=binomial(link='logit'),data=Soil_characteristics)
summary(GLMER_soil)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Location ~ scale(Agg_stability) * scale(C_percent) * scale(N_percent) +
## (1 | Population)
## Data: Soil_characteristics
##
## AIC BIC logLik deviance df.resid
## 36.4 45.3 -9.2 18.4 11
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4454 -0.7165 -0.1067 0.6076 3.1160
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 6.629e-17 8.142e-09
## Number of obs: 20, groups: Population, 10
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 1.895 1.494
## scale(Agg_stability) -3.634 2.731
## scale(C_percent) 10.664 7.517
## scale(N_percent) -10.010 7.247
## scale(Agg_stability):scale(C_percent) -15.686 10.695
## scale(Agg_stability):scale(N_percent) 15.667 10.403
## scale(C_percent):scale(N_percent) -1.937 1.343
## scale(Agg_stability):scale(C_percent):scale(N_percent) 3.464 2.511
## z value Pr(>|z|)
## (Intercept) 1.269 0.204
## scale(Agg_stability) -1.331 0.183
## scale(C_percent) 1.419 0.156
## scale(N_percent) -1.381 0.167
## scale(Agg_stability):scale(C_percent) -1.467 0.142
## scale(Agg_stability):scale(N_percent) 1.506 0.132
## scale(C_percent):scale(N_percent) -1.442 0.149
## scale(Agg_stability):scale(C_percent):scale(N_percent) 1.380 0.168
##
## Correlation of Fixed Effects:
## (Intr) sc(A_) sc(C_) sc(N_) sc(A_):(C_) s(A_):(N s(C_):
## scl(Agg_st) -0.740
## scl(C_prcn) 0.617 -0.441
## scl(N_prcn) -0.522 0.384 -0.986
## sc(A_):(C_) -0.608 0.791 -0.637 0.604
## sc(A_):(N_) 0.558 -0.700 0.635 -0.602 -0.981
## sc(C_):(N_) -0.829 0.730 -0.454 0.350 0.695 -0.678
## s(A_):(C_): 0.718 -0.866 0.656 -0.629 -0.774 0.684 -0.635options(na.action = "na.fail")
dredge(GLMER_soil) ## Fixed term is "(Intercept)"# A model with nothing is best...Histology data Analysis
Histology Data exploration
Calculate plant-level Mycchorizae, pathogens and herbivory:
# Fix column characteristics for plotting and analysis
fdata$indiv2<-as.factor(fdata$indiv2)
fdata$Mycorrhiza<-as.numeric(paste(fdata$Mycorrhiza))
fdata$Lesion<-as.numeric(paste(fdata$Lesion))
fdata$Herbivory<-as.numeric(paste(fdata$Herbivory))
pct<-function(x){
sum(x)/100
}
# Calculate plant-level averages
fdat_ind<-aggregate(fdata[,c("Mycorrhiza","Lesion","Herbivory")],by=list(indiv2=fdata$indiv2,pop=fdata$pop,loc=fdata$loc,species=fdata$species), FUN=pct)Simple linear model visualization:
## Mycorrhiza model (NOTE: No loc*species because not all species present in all populations)
base.Myc<-lmer(Mycorrhiza~species+loc+species:loc+(1|pop)+(1|pop:loc),data=fdat_ind)
noint.Myc<-lmer(Mycorrhiza~species+loc+(1|pop)+(1|pop:loc),data=fdat_ind)
anova(base.Myc,noint.Myc)## refitting model(s) with ML (instead of REML)base.Myc ## Inspect coefficients## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Mycorrhiza ~ species + loc + species:loc + (1 | pop) + (1 | pop:loc)
## Data: fdat_ind
## REML criterion at convergence: -12.5748
## Random effects:
## Groups Name Std.Dev.
## pop:loc (Intercept) 0.06622
## pop (Intercept) 0.00000
## Residual 0.21042
## Number of obs: 225, groups: pop:loc, 22; pop, 11
## Fixed Effects:
## (Intercept) speciesCL speciesGA speciesGR
## 0.38985 -0.05067 0.03105 0.01338
## speciesMR speciesSC locO speciesCL:locO
## -0.02840 0.28780 0.15278 0.13277
## speciesGA:locO speciesGR:locO speciesMR:locO speciesSC:locO
## -0.10069 0.23077 0.28893 -0.04707ggplot(fdat_ind,aes(species,Mycorrhiza,colour=loc))+geom_boxplot()+theme_simple()
ggplot(fdat_ind,aes(pop,Mycorrhiza,colour=loc))+geom_boxplot()+theme_simple()
Mycorrhia model results
- Mycorrhiza colonization is significantly higher outside patches of A. petiolata than inside 2a. Species vary in Mycorrhiza colonization, with highest colonization in SC and others relatively similar. 2b. Effect of A. petiolata on Mycorrhiza colonization differs among species, with more pronounced effects for GR and MR, and will little effect for GA
- Effect of A. petiolata on Mycorrhiza colonization varies by population, but not as strongly as among species
## Mycorrhiza model (NOTE: No loc*species because not all species present in all populations)
base.Les<-lmer(Lesion~Mycorrhiza*species*loc+(1|pop/loc),data=fdat_ind)
norint.Les<-lmer(Lesion~Mycorrhiza*species*loc+(1|pop),data=fdat_ind)
anova(base.Les,norint.Les)## refitting model(s) with ML (instead of REML)norint.Les # inspect coefficients## Linear mixed model fit by REML ['lmerMod']
## Formula: Lesion ~ Mycorrhiza * species * loc + (1 | pop)
## Data: fdat_ind
## REML criterion at convergence: -70.9219
## Random effects:
## Groups Name Std.Dev.
## pop (Intercept) 0.04934
## Residual 0.18543
## Number of obs: 225, groups: pop, 11
## Fixed Effects:
## (Intercept) Mycorrhiza
## 0.686711 0.053896
## speciesCL speciesGA
## -0.188794 0.146342
## speciesGR speciesMR
## 0.090499 -0.310235
## speciesSC locO
## -0.037506 0.016259
## Mycorrhiza:speciesCL Mycorrhiza:speciesGA
## 0.216958 -0.702730
## Mycorrhiza:speciesGR Mycorrhiza:speciesMR
## -0.020628 0.559261
## Mycorrhiza:speciesSC Mycorrhiza:locO
## 0.136864 -0.056951
## speciesCL:locO speciesGA:locO
## -0.225004 -0.389157
## speciesGR:locO speciesMR:locO
## -0.565486 -0.129568
## speciesSC:locO Mycorrhiza:speciesCL:locO
## -0.007965 0.098911
## Mycorrhiza:speciesGA:locO Mycorrhiza:speciesGR:locO
## 0.784046 0.433494
## Mycorrhiza:speciesMR:locO Mycorrhiza:speciesSC:locO
## -0.172551 -0.051787Les2<-lmer(Lesion~Mycorrhiza*species*loc-Mycorrhiza:species:loc+(1|pop),data=fdat_ind)
anova(norint.Les,Les2)## refitting model(s) with ML (instead of REML)Les3a<-lmer(Lesion~Mycorrhiza*species+Mycorrhiza*loc+(1|pop),data=fdat_ind)
Les3b<-lmer(Lesion~Mycorrhiza*species+species*loc+(1|pop),data=fdat_ind)
Les3c<-lmer(Lesion~species*loc+Mycorrhiza*loc+(1|pop),data=fdat_ind)
anova(norint.Les,Les3a) # TEST species*loc ## refitting model(s) with ML (instead of REML)anova(norint.Les,Les3b) # TEST Myc*loc ## refitting model(s) with ML (instead of REML)anova(norint.Les,Les3c) # TEST Myc*species## refitting model(s) with ML (instead of REML)ggplot(fdat_ind,aes(Mycorrhiza,Lesion))+facet_wrap(~species)+
geom_point(aes(col=loc))+theme_simple()+geom_smooth(method="lm")
** Lesions model results**
- Root lesions increase with mycorrhiza colonization, contrary to predictions
- Effect of Mycorrhiza on Lesions depends strongly on species (Myc*species interactions)
Hypothesis 1: Lesion ~ Mycorriza
GIVEN THE NEW ANALYSIS ABOVE, THIS ISN’T WORTHWILE BECAUSE THERE IS NO SIGNIFICANT Myc:species:loc
# calculating the r
indivCODE <- unique(fdata$indiv2)
dfcorr <-data.frame(indivCODE=c(0),r=c(0))
for (i in 1:length(indivCODE)){
dat1<-NULL
dat1<-filter(fdata, fdata$indiv2==indivCODE[i])
dfcorr[i,1] <- indivCODE[i]
dfcorr[i,2] <- cor(as.numeric(dat1$Mycorrhiza),as.numeric(dat1$Lesion))
}## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zerodfcorr %>% drop_na()->dfcorrCo-occurence data
#Normal?
car::qqp(dfcorr$r, "norm")
shapiro.test(dfcorr$r)
hist(dfcorr$r)
# Cannot be normal because it is only between -1 and 1 and our mean is at 0.
descdist(dfcorr$r, discrete = FALSE, boot = 500)
# weibull when transformed to be positive?
dfcorr$r.t <- dfcorr$r+1
fit.weibull <- fitdist(dfcorr$r.t, "weibull")
car::qqp(dfcorr$r.t, "weibull", shape = fit.weibull$estimate[[1]])
dfcorr$pop <- gsub("^([0-9]+)[I,O][A-Z]+.*","\\1",dfcorr$indivCODE)
dfcorr$location <- gsub("^[0-9]+([I,O])[A-Z]+.*","\\1",dfcorr$indivCODE)
dfcorr$species <- sub("^[0-9]+[I,O]([A-Z]+).*","\\1",dfcorr$indivCODE)
dfcorr$location <- as.factor(dfcorr$location)
dfcorr$species <- as.factor(dfcorr$species)
dfcorr$pop <- as.factor(dfcorr$pop)
meanR<-mean(dfcorr$r)
SEMR<-sd(dfcorr$r)/sqrt(length(dfcorr$r))
UCI<- meanR+SEMR*1.96
LCI<- meanR-SEMR*1.96
mean(filter(dfcorr, dfcorr$location=="I")$r)
mean(filter(dfcorr, dfcorr$location=="O")$r)
baysian_corrGLMM <- brm(r.t ~ location * species + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM2 <- brm(r.t ~ location + species + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM3 <- brm(r.t ~ location + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM_NULL <- brm(r.t ~ 1 + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
waic.BCG_1<-waic(baysian_corrGLMM)
waic.BCG_2<-waic(baysian_corrGLMM2)
waic.BCG_3<-waic(baysian_corrGLMM3)
waic.BCG_NULL<-waic(baysian_corrGLMM_NULL)
compare(waic.BCG_1,waic.BCG_2) #positive means second is better
compare(waic.BCG_2,waic.BCG_3) #negative means first is better
compare(waic.BCG_2,waic.BCG_NULL) #negative means first is better
summary(baysian_corrGLMM2)
bayes_R2(baysian_corrGLMM2)
ggplot(dfcorr, aes(y=r,x=location))+geom_boxplot()+theme_simple()
ggplot(dfcorr, aes(y=r,x=species))+geom_boxplot()+theme_simple()Totals data
# model
GLMMmyc_lesion <- glmer(binlesion ~ None.Myc * species + (1 | pop) , data = ndatas,family = binomial(link = "logit"))
summary(GLMMmyc_lesion)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: binlesion ~ None.Myc * species + (1 | pop)
## Data: ndatas
##
## AIC BIC logLik deviance df.resid
## 5388.9 5433.4 -2681.5 5362.9 213
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9159 -3.2800 -0.7426 2.8805 11.6817
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop (Intercept) 0.07038 0.2653
## Number of obs: 226, groups: pop, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.750276 0.090719 -8.270 < 2e-16 ***
## None.Myc -0.017781 0.038118 -0.466 0.640866
## speciesCL 0.415585 0.056183 7.397 1.39e-13 ***
## speciesGA 0.444280 0.053486 8.306 < 2e-16 ***
## speciesGR 0.266576 0.073678 3.618 0.000297 ***
## speciesMR 0.146355 0.057804 2.532 0.011344 *
## speciesSC 0.007849 0.051574 0.152 0.879040
## None.Myc:speciesCL 0.156587 0.048414 3.234 0.001219 **
## None.Myc:speciesGA 0.191932 0.050731 3.783 0.000155 ***
## None.Myc:speciesGR 0.338206 0.069696 4.853 1.22e-06 ***
## None.Myc:speciesMR 0.063131 0.053541 1.179 0.238360
## None.Myc:speciesSC 0.210590 0.051242 4.110 3.96e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Nn.Myc spcsCL spcsGA spcsGR spcsMR spcsSC N.M:CL N.M:GA
## None.Myc -0.134
## speciesCL -0.369 0.213
## speciesGA -0.385 0.223 0.649
## speciesGR -0.279 0.160 0.466 0.500
## speciesMR -0.340 0.207 0.572 0.614 0.438
## speciesSC -0.381 0.233 0.627 0.660 0.477 0.588
## Nn.Myc:spCL 0.098 -0.786 -0.235 -0.155 -0.110 -0.155 -0.171
## Nn.Myc:spGA 0.094 -0.750 -0.141 -0.235 -0.116 -0.146 -0.169 0.593
## Nn.Myc:spGR 0.086 -0.551 -0.142 -0.151 -0.111 -0.128 -0.145 0.431 0.406
## Nn.Myc:spMR 0.090 -0.711 -0.144 -0.148 -0.092 -0.111 -0.159 0.560 0.534
## Nn.Myc:spSC 0.116 -0.750 -0.191 -0.195 -0.145 -0.169 -0.045 0.589 0.555
## N.M:GR N.M:MR
## None.Myc
## speciesCL
## speciesGA
## speciesGR
## speciesMR
## speciesSC
## Nn.Myc:spCL
## Nn.Myc:spGA
## Nn.Myc:spGR
## Nn.Myc:spMR 0.394
## Nn.Myc:spSC 0.429 0.527options(na.action="na.fail")
dredge(GLMMmyc_lesion)## Fixed term is "(Intercept)"r.squaredGLMM(GLMMmyc_lesion)## The result is correct only if all data used by the model has not changed since model was fitted.## R2m R2c
## 0.02739802 0.04008648# will need to change this ASAP
ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path))+
stat_summary(fun.y=mean, geom="point")+
geom_smooth(method = "lm")+
geom_point()+
theme_simple()+
facet_wrap("species")+
labs(x = "Mycorrhizal colonization", y="Lesions")
Hypothesis 2: Mycorrhiza & Lesion ~ invasion
Does location impact the effect of myc on lesion?
GLMMmyc_lesion_loc <- glmer(binlesion ~ None.Myc * location * species + (1 | pop), data = ndatas, family = binomial(link = "logit"))
summary(GLMMmyc_lesion_loc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: binlesion ~ None.Myc * location * species + (1 | pop)
## Data: ndatas
##
## AIC BIC logLik deviance df.resid
## 5185.5 5271.0 -2567.7 5135.5 201
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.6366 -3.2043 -0.3462 2.6760 13.6450
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop (Intercept) 0.0698 0.2642
## Number of obs: 226, groups: pop, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.725383 0.100890 -7.190 6.49e-13 ***
## None.Myc -0.062764 0.054987 -1.141 0.253688
## locationO -0.068081 0.081316 -0.837 0.402458
## speciesCL 0.200476 0.081638 2.456 0.014063 *
## speciesGA 0.296660 0.077883 3.809 0.000139 ***
## speciesGR -1.181296 0.193649 -6.100 1.06e-09 ***
## speciesMR 0.277361 0.084158 3.296 0.000982 ***
## speciesSC -0.021562 0.073106 -0.295 0.768042
## None.Myc:locationO 0.093031 0.080266 1.159 0.246440
## None.Myc:speciesCL 0.324299 0.071788 4.517 6.26e-06 ***
## None.Myc:speciesGA 0.189403 0.072737 2.604 0.009216 **
## None.Myc:speciesGR -0.300771 0.187873 -1.601 0.109393
## None.Myc:speciesMR 0.255496 0.077895 3.280 0.001038 **
## None.Myc:speciesSC 0.101743 0.070113 1.451 0.146743
## locationO:speciesCL 0.398803 0.106396 3.748 0.000178 ***
## locationO:speciesGA 0.301929 0.101835 2.965 0.003028 **
## locationO:speciesGR 1.889287 0.212942 8.872 < 2e-16 ***
## locationO:speciesMR -0.260917 0.111010 -2.350 0.018754 *
## locationO:speciesSC 0.143608 0.103149 1.392 0.163850
## None.Myc:locationO:speciesCL -0.282839 0.100392 -2.817 0.004842 **
## None.Myc:locationO:speciesGA -0.003347 0.105512 -0.032 0.974693
## None.Myc:locationO:speciesGR 0.463051 0.207868 2.228 0.025906 *
## None.Myc:locationO:speciesMR -0.384191 0.109499 -3.509 0.000450 ***
## None.Myc:locationO:speciesSC 0.233745 0.102953 2.270 0.023183 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itemm<-emmeans(GLMMmyc_lesion_loc, list(~None.Myc| species | location, ~None.Myc| location), type = "response" )## NOTE: Results may be misleading due to involvement in interactionssummary(emm)## $`emmeans of None.Myc, location | species`
## species = AH:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.3262087 0.02217532 NA 0.28432176 0.3710671
## 1.058488e-16 O 0.3114255 0.02097601 NA 0.27186461 0.3539447
##
## species = CL:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.3717055 0.02222786 NA 0.32927810 0.4162074
## 1.058488e-16 O 0.4516058 0.02288001 NA 0.40727139 0.4967219
##
## species = GA:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.3944313 0.02175155 NA 0.35269620 0.4377651
## 1.058488e-16 O 0.4514349 0.02239380 NA 0.40803107 0.4955905
##
## species = GR:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.1293544 0.02255808 NA 0.09118441 0.1803321
## 1.058488e-16 O 0.4786449 0.02657122 NA 0.42698389 0.5307668
##
## species = MR:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.3898311 0.02331380 NA 0.34521822 0.4363672
## 1.058488e-16 O 0.3149625 0.02075782 NA 0.27576180 0.3569890
##
## species = SC:
## None.Myc location prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 I 0.3214874 0.01942629 NA 0.28465376 0.3606841
## 1.058488e-16 O 0.3381795 0.02091737 NA 0.29847025 0.3803078
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $`emmeans of None.Myc | location`
## location = I:
## None.Myc prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 0.3107155 0.01873775 NA 0.2752289 0.3485770
##
## location = O:
## None.Myc prob SE df asymp.LCL asymp.UCL
## 1.058488e-16 0.3887098 0.01966447 NA 0.3509329 0.4278726
##
## Results are averaged over the levels of: species
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scaleggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path))+
geom_smooth(method = "lm")+
geom_point(aes(color=location))+
theme_simple()+
facet_wrap("species")+
labs(x = "Mycorrhizal colonization", y="Lesions")
pca <- prcomp(ndata[,c(5,9)])
location <- ndata[,13]
g <- ggbiplot(pca, obs.scale = 1, var.scale = 1,
groups = location, ellipse = TRUE,
circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g)
Does location impact lesion alone?
glmerLESION_disp <- glmer (binlesion ~ location * species + (1|pop) + (1|indiv), data = ndatas, family = binomial(link="logit"))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00135513 (tol =
## 0.001, component 1)summary(glmerLESION_disp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: binlesion ~ location * species + (1 | pop) + (1 | indiv)
## Data: ndatas
##
## AIC BIC logLik deviance df.resid
## 2011.7 2059.6 -991.8 1983.7 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.53231 -0.13210 -0.00971 0.11216 1.14859
##
## Random effects:
## Groups Name Variance Std.Dev.
## indiv (Intercept) 1.11959 1.0581
## pop (Intercept) 0.06426 0.2535
## Number of obs: 226, groups: indiv, 226; pop, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.85482 0.30710 -2.784 0.00538 **
## locationO -0.10392 0.39475 -0.263 0.79235
## speciesCL 0.26989 0.39300 0.687 0.49225
## speciesGA 0.38852 0.37584 1.034 0.30127
## speciesGR -0.88023 0.62967 -1.398 0.16214
## speciesMR 0.38516 0.41882 0.920 0.35776
## speciesSC -0.09661 0.35549 -0.272 0.78581
## locationO:speciesCL 0.45477 0.52759 0.862 0.38870
## locationO:speciesGA 0.47505 0.50259 0.945 0.34456
## locationO:speciesGR 1.77188 0.77503 2.286 0.02224 *
## locationO:speciesMR -0.28686 0.55694 -0.515 0.60651
## locationO:speciesSC -0.02068 0.48262 -0.043 0.96582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) loctnO spcsCL spcsGA spcsGR spcsMR spcsSC lcO:CL lcO:GA
## locationO -0.703
## speciesCL -0.741 0.551
## speciesGA -0.778 0.579 0.616
## speciesGR -0.455 0.347 0.354 0.376
## speciesMR -0.682 0.522 0.539 0.570 0.340
## speciesSC -0.807 0.608 0.637 0.669 0.393 0.588
## lctnO:spcCL 0.524 -0.747 -0.716 -0.431 -0.253 -0.390 -0.453
## lctnO:spcGA 0.566 -0.787 -0.446 -0.727 -0.272 -0.411 -0.488 0.587
## lctnO:spcGR 0.354 -0.513 -0.273 -0.288 -0.807 -0.269 -0.307 0.378 0.400
## lctnO:spcMR 0.502 -0.709 -0.396 -0.411 -0.250 -0.732 -0.434 0.530 0.558
## lctnO:spcSC 0.575 -0.818 -0.450 -0.473 -0.283 -0.428 -0.721 0.612 0.643
## lcO:GR lcO:MR
## locationO
## speciesCL
## speciesGA
## speciesGR
## speciesMR
## speciesSC
## lctnO:spcCL
## lctnO:spcGA
## lctnO:spcGR
## lctnO:spcMR 0.366
## lctnO:spcSC 0.419 0.580
## convergence code: 0
## Model failed to converge with max|grad| = 0.00135513 (tol = 0.001, component 1)options(na.action="na.fail")
dredge(glmerLESION_disp)## Fixed term is "(Intercept)"## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00135513 (tol =
## 0.001, component 1)r.squaredGLMM(glmerLESION_disp)## R2m R2c
## 0.03076813 0.04469034ggplot(aes(x = location, y = 100-None.Path, color=location), data=ndata)+
geom_boxplot()+
geom_jitter(width=0.1)+
facet_wrap("species")+
theme_simple()
Does location impact mycorrhiza alone?
glmerMYC_disp <- glmer (binmycorr ~ location * species + (1|pop) + (1|indiv), data = ndata, family = binomial(link="logit"))
summary(glmerMYC_disp)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: binmycorr ~ location * species + (1 | pop) + (1 | indiv)
## Data: ndata
##
## AIC BIC logLik deviance df.resid
## 2073.3 2121.2 -1022.7 2045.3 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.35863 -0.12130 0.01523 0.09795 1.40036
##
## Random effects:
## Groups Name Variance Std.Dev.
## indiv (Intercept) 2.062e+00 1.4361234
## pop (Intercept) 1.044e-08 0.0001022
## Number of obs: 226, groups: indiv, 226; pop, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.121220 0.390387 0.310 0.7562
## locationO -0.008745 0.526879 -0.017 0.9868
## speciesCL -0.200766 0.515746 -0.389 0.6971
## speciesGA -0.375379 0.491194 -0.764 0.4447
## speciesGR -1.324367 0.829848 -1.596 0.1105
## speciesMR -0.505583 0.551771 -0.916 0.3595
## speciesSC -0.954255 0.471189 -2.025 0.0428 *
## locationO:speciesCL 0.053276 0.706272 0.075 0.9399
## locationO:speciesGA 0.233201 0.669292 0.348 0.7275
## locationO:speciesGR 0.990565 1.025133 0.966 0.3339
## locationO:speciesMR -0.245112 0.745205 -0.329 0.7422
## locationO:speciesSC -0.457001 0.645298 -0.708 0.4788
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) loctnO spcsCL spcsGA spcsGR spcsMR spcsSC lcO:CL lcO:GA
## locationO -0.741
## speciesCL -0.757 0.561
## speciesGA -0.795 0.589 0.602
## speciesGR -0.470 0.349 0.356 0.374
## speciesMR -0.708 0.524 0.536 0.562 0.333
## speciesSC -0.829 0.614 0.627 0.659 0.390 0.586
## lctnO:spcCL 0.553 -0.746 -0.730 -0.439 -0.260 -0.391 -0.458
## lctnO:spcGA 0.583 -0.787 -0.442 -0.734 -0.274 -0.413 -0.483 0.587
## lctnO:spcGR 0.381 -0.514 -0.288 -0.303 -0.809 -0.269 -0.316 0.383 0.405
## lctnO:spcMR 0.524 -0.707 -0.396 -0.416 -0.246 -0.740 -0.434 0.527 0.557
## lctnO:spcSC 0.605 -0.816 -0.458 -0.481 -0.285 -0.428 -0.730 0.609 0.643
## lcO:GR lcO:MR
## locationO
## speciesCL
## speciesGA
## speciesGR
## speciesMR
## speciesSC
## lctnO:spcCL
## lctnO:spcGA
## lctnO:spcGR
## lctnO:spcMR 0.363
## lctnO:spcSC 0.420 0.577options(na.action="na.fail")
dredge(glmerMYC_disp)## Fixed term is "(Intercept)"## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0029301 (tol =
## 0.001, component 1)r.squaredGLMM(glmerMYC_disp)## R2m R2c
## 0.0410656 0.0410656ggplot(aes(x = location, y = 100-None.Myc, color=location), data=ndata)+
geom_boxplot()+
geom_jitter(width=0.1)+
facet_wrap("species")+
theme_simple()
Hypothesis 3 + 4: Lesion ~ Mycorriza * invasion
ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path, color=location))+
geom_smooth(method = "lm")+
geom_point()+
theme_simple()+
facet_wrap("species")+
labs(x = "Mycorrhizal colonization", y="Lesions")
ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path, color=location))+
geom_smooth(method = "lm")+
geom_point()+
theme_simple()+
labs(x = "Mycorrhizal colonization", y="Lesions")
Bring things together
Soil_characteristics$Location<-as.character(Soil_characteristics$Location)
for(i in 1:length(Soil_characteristics$Location))if(Soil_characteristics[i,2]=="in"){
Soil_characteristics[i,2]<-"I"
}else{
Soil_characteristics[i,2]<-"O"
}
Soil_characteristics$Location<-as.factor(Soil_characteristics$Location)
names(Soil_characteristics)[1:2]<-c("pop","location")
fulldata <- merge(Soil_characteristics, ndata, by=c("pop","location"))
# note i used the SCALED dataframe for the total histology
fulldata <- merge(fulldata, Root_size, by=c("indiv","location","pop","species"))
for(i in 1:length(fulldata$indiv)){
fulldata$mean_rootsize[i]<-(sum(fulldata[i,20:29])/length(fulldata[,20:29]))
}
fulldata<-fulldata[,-c(20:29)]
binlesionF<-cbind(fulldata$None.Path,100-fulldata$None.Path)
binmycorrF<-cbind(fulldata$None.Myc,100-fulldata$None.Myc)
binherbivoryF<-cbind(fulldata$None.Herb,100-fulldata$None.Herb)
fulldata$pop<-as.factor(fulldata$pop)
# model myc_lesion
mod1 <- glmmadmb(binlesionF ~ scale(None.Myc) * scale(mean_rootsize) * species + scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) + (1 | pop), data = fulldata,family = "binomial")
summary(mod1)##
## Call:
## glmmadmb(formula = binlesionF ~ scale(None.Myc) * scale(mean_rootsize) *
## species + scale(pH) + scale(Agg_stability) + scale(N_percent) +
## scale(C_percent) + (1 | pop), data = fulldata, family = "binomial")
##
## AIC: 3914.8
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -0.9247 0.1360 -6.80
## scale(None.Myc) -0.0167 0.0601 -0.28
## scale(mean_rootsize) 0.1370 0.0727 1.88
## speciesCL 0.4485 0.0730 6.15
## speciesGA 1.7767 0.1557 11.41
## speciesGR 0.3890 0.1028 3.78
## speciesMR 0.3305 0.0742 4.45
## speciesSC 0.4228 0.0721 5.87
## scale(pH) -0.1158 0.0262 -4.43
## scale(Agg_stability) -0.0166 0.0322 -0.51
## scale(N_percent) -0.6563 0.1253 -5.24
## scale(C_percent) 0.6935 0.1489 4.66
## scale(None.Myc):scale(mean_rootsize) 0.1930 0.0619 3.12
## scale(None.Myc):speciesCL 0.0459 0.0712 0.64
## scale(None.Myc):speciesGA 0.3262 0.1761 1.85
## scale(None.Myc):speciesGR 0.0850 0.1080 0.79
## scale(None.Myc):speciesMR 0.1156 0.0755 1.53
## scale(None.Myc):speciesSC 0.3634 0.0793 4.58
## scale(mean_rootsize):speciesCL 0.0602 0.0806 0.75
## scale(mean_rootsize):speciesGA 1.0349 0.1636 6.33
## scale(mean_rootsize):speciesGR -0.3641 0.1051 -3.46
## scale(mean_rootsize):speciesMR -0.0544 0.0882 -0.62
## scale(mean_rootsize):speciesSC 0.0131 0.1275 0.10
## scale(None.Myc):scale(mean_rootsize):speciesCL -0.0199 0.0656 -0.30
## scale(None.Myc):scale(mean_rootsize):speciesGA -0.1636 0.1831 -0.89
## scale(None.Myc):scale(mean_rootsize):speciesGR -0.4582 0.1258 -3.64
## scale(None.Myc):scale(mean_rootsize):speciesMR -0.2581 0.0762 -3.39
## scale(None.Myc):scale(mean_rootsize):speciesSC -0.1735 0.1278 -1.36
## Pr(>|z|)
## (Intercept) 1.1e-11 ***
## scale(None.Myc) 0.78093
## scale(mean_rootsize) 0.05950 .
## speciesCL 7.9e-10 ***
## speciesGA < 2e-16 ***
## speciesGR 0.00016 ***
## speciesMR 8.4e-06 ***
## speciesSC 4.5e-09 ***
## scale(pH) 9.6e-06 ***
## scale(Agg_stability) 0.60758
## scale(N_percent) 1.6e-07 ***
## scale(C_percent) 3.2e-06 ***
## scale(None.Myc):scale(mean_rootsize) 0.00181 **
## scale(None.Myc):speciesCL 0.51893
## scale(None.Myc):speciesGA 0.06396 .
## scale(None.Myc):speciesGR 0.43096
## scale(None.Myc):speciesMR 0.12571
## scale(None.Myc):speciesSC 4.5e-06 ***
## scale(mean_rootsize):speciesCL 0.45510
## scale(mean_rootsize):speciesGA 2.5e-10 ***
## scale(mean_rootsize):speciesGR 0.00053 ***
## scale(mean_rootsize):speciesMR 0.53735
## scale(mean_rootsize):speciesSC 0.91845
## scale(None.Myc):scale(mean_rootsize):speciesCL 0.76130
## scale(None.Myc):scale(mean_rootsize):speciesGA 0.37155
## scale(None.Myc):scale(mean_rootsize):speciesGR 0.00027 ***
## scale(None.Myc):scale(mean_rootsize):speciesMR 0.00071 ***
## scale(None.Myc):scale(mean_rootsize):speciesSC 0.17451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of observations: total=189, pop=9
## Random effect variance(s):## Warning in .local(x, sigma, ...): 'sigma' and 'rdig' arguments are present
## for compatibility only: ignored## Group=pop
## Variance StdDev
## (Intercept) 0.1348 0.3671
##
##
## Log-likelihood: -1928.4mod1.dredge <- dredge(mod1)## Fixed term is "(Intercept)"mod1.dredgegood <- subset(mod1.dredge,mod1.dredge$delta<=7)
# model myc
mod2 <- glmmadmb(binmycorrF ~ scale(mean_rootsize) * species + scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) + (1 | pop) , data = fulldata,family = "binomial")
summary(mod2)##
## Call:
## glmmadmb(formula = binmycorrF ~ scale(mean_rootsize) * species +
## scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) +
## (1 | pop), data = fulldata, family = "binomial")
##
## AIC: 6440.8
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.00125 0.10116 -0.01 0.99014
## scale(mean_rootsize) 0.12660 0.06462 1.96 0.05008 .
## speciesCL -0.04954 0.06748 -0.73 0.46286
## speciesGA -0.03748 0.13458 -0.28 0.78065
## speciesGR -0.40656 0.08369 -4.86 1.2e-06 ***
## speciesMR -0.25895 0.06867 -3.77 0.00016 ***
## speciesSC -0.82100 0.06393 -12.84 < 2e-16 ***
## scale(pH) 0.08393 0.02423 3.46 0.00053 ***
## scale(Agg_stability) 0.18425 0.03032 6.08 1.2e-09 ***
## scale(N_percent) 0.14058 0.10894 1.29 0.19691
## scale(C_percent) -0.14612 0.13142 -1.11 0.26619
## scale(mean_rootsize):speciesCL -0.17833 0.07173 -2.49 0.01292 *
## scale(mean_rootsize):speciesGA -0.17024 0.14277 -1.19 0.23309
## scale(mean_rootsize):speciesGR -0.26173 0.08015 -3.27 0.00109 **
## scale(mean_rootsize):speciesMR -0.12646 0.07473 -1.69 0.09062 .
## scale(mean_rootsize):speciesSC -1.26599 0.12042 -10.51 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of observations: total=189, pop=9
## Random effect variance(s):## Warning in .local(x, sigma, ...): 'sigma' and 'rdig' arguments are present
## for compatibility only: ignored## Group=pop
## Variance StdDev
## (Intercept) 0.0652 0.2554
##
##
## Log-likelihood: -3203.42mod2.dredge <- dredge(mod2)## Fixed term is "(Intercept)"mod2.dredgegood <- subset(mod2.dredge,mod2.dredge$delta<=7)